1 DATA

## Longitudinal data
data_sum <- loadlongitudinaldata(dataset = "DATA_Adults_G1G29.csv", rm_generation1 = 1,rm_generation2 = 7,rm_generation3 = 29)

## Phenotyping steps
data_G0 <- loadfitnessdata(dataset = "Selection_Phenotypage_G0_G7_G8.csv", generation = "G1")
data_G7 <- loadfitnessdata(dataset = "Selection_Phenotypage_G0_G7_G8.csv", generation = "G7")
data_G29 <- loadfitnessdata(dataset = 
                            "PERFORMANCE_Comptage_adultes_G13G14G15G16G17G18G19G20G21G22G23G24G25G26G27G28G29.csv",
                            generation = "29")
head(data_sum)
##   Line Fruit_s Generation         Phase  N Nb_adults       sd    fitness
## 1  CE1  Cherry          2 first_prepool 20   9.15000 6.123939 -0.7819784
## 2  CE1  Cherry          3 first_prepool 10  15.30000 7.631077 -0.2678794
## 3  CE1  Cherry          4 first_prepool  8  15.00000 7.782765 -0.2876821
## 4  CE1  Cherry          5 first_prepool  6  14.50000 6.284903 -0.3215836
## 5  CE2  Cherry          2 first_prepool 20   7.75000 7.758696 -0.9480394
## 6  CE2  Cherry          3 first_prepool  7  14.42857 8.303757 -0.3265219
##   se_fitness
## 1  0.1496562
## 2  0.1577228
## 3  0.1834415
## 4  0.1769518
## 5  0.2238577
## 6  0.2175215
head(data_G0)
##     Treatment Line Fruit_s Nb_eggs Nb_adults SA Emergence_rate
## 993    Cherry  CE1      GF      76         6  0     0.07894737
## 994    Cherry  CE1      GF      89        17  0     0.19101124
## 995    Cherry  CE1      GF      57        12  0     0.21052632
## 996    Cherry  CE1      GF     172        24  0     0.13953488
## 997    Cherry  CE1      GF     173        33  0     0.19075145
## 998    Cherry  CE1      GF      91        18  0     0.19780220
head(data_G7)
##    Treatment Line    Fruit_s Nb_eggs Nb_adults SA Emergence_rate
## 3 Strawberry  CR4  Cranberry     152        68  0      0.4473684
## 4  Cranberry  CR4  Cranberry     246        25  1      0.1016260
## 5     Cherry  CR4  Cranberry     238        29  0      0.1218487
## 6     Cherry  CR4  Cranberry     166        23  0      0.1385542
## 8  Cranberry  FR3 Strawberry     204         5  0      0.0245098
## 9 Strawberry  FR3 Strawberry     124        45  1      0.3629032
head(data_G29)
##       Treatment Line Fruit_s Nb_eggs Nb_adults SA Emergence_rate
## 5392 Strawberry  CEA  Cherry     196        16  0     0.08163265
## 5393 Strawberry  CEA  Cherry     192        30  0     0.15625000
## 5394 Strawberry  CEA  Cherry     160        17  0     0.10625000
## 5395 Strawberry  CEA  Cherry     106         9  0     0.08490566
## 5396 Strawberry  CEA  Cherry     119        14  0     0.11764706
## 5397 Strawberry  CEA  Cherry     204        24  0     0.11764706
dim(data_G29)
## [1] 990   7
#Compute log change
## !! sd_logchange are standard errors, not standard deviations
data_logchange <- rbind(computelogchangenew(fitness_dataset_ancestral = data_G0, fitness_dataset_evolved = data_G7, trait = "Nb_adults", gen=7), computelogchangenew(fitness_dataset_ancestral = data_G0, fitness_dataset_evolved = data_G29, trait = "Nb_adults", gen=29))

data_logchangeNb_eggs <- rbind(computelogchangenew(fitness_dataset_ancestral = data_G0, fitness_dataset_evolved = data_G7, trait = "Nb_eggs", gen=7), computelogchangenew(fitness_dataset_ancestral = data_G0, fitness_dataset_evolved = data_G29, trait = "Nb_eggs", gen=29))

data_logchangeEmergence_rate <- rbind(computelogchangenew(fitness_dataset_ancestral = data_G0, fitness_dataset_evolved = data_G7, trait = "Emergence_rate", gen=7), computelogchangenew(fitness_dataset_ancestral = data_G0, fitness_dataset_evolved = data_G29, trait = "Emergence_rate", gen=29))

## Combine datasets
data_logchangeNb_eggsEmergence_rate <- data.frame(data_logchange, Nb_eggslogchange=data_logchangeNb_eggs$logchange, sdNb_eggslogchange=data_logchangeNb_eggs$sd_logchange, Emergence_ratelogchange=data_logchangeEmergence_rate$logchange, sdEmergence_ratelogchange=data_logchangeEmergence_rate$sd_logchange)

#Lines: CE2 and CR2 do not pass Geary's test (the threshold of a standardized mean greater than 3)
data_logchange$Line_type  <- ifelse(data_logchange$Line == "CR2"| data_logchange$Line == "CE2", "solid","dashed")
data_logchangeNb_eggs$Line_type  <- ifelse(data_logchange$Line == "CR2"| data_logchange$Line == "CE2", "solid","dashed")
data_logchangeEmergence_rate$Line_type  <- ifelse(data_logchange$Line == "CR2"| data_logchange$Line == "CE2", "solid","dashed")


#Formatting for testing correlation
TEMP_dataG7_CheCran <- formattinglogchange(data_logchange, "7", "Cherry", "Cranberry")
TEMP_dataG7_CranStraw <- formattinglogchange(data_logchange, "7", "Cranberry", "Strawberry")
TEMP_dataG7_StrawChe <- formattinglogchange(data_logchange, "7", "Strawberry", "Cherry")

TEMP_dataG29_CheCran <- formattinglogchange(data_logchange, "29", "Cherry", "Cranberry")
TEMP_dataG29_CranStraw <- formattinglogchange(data_logchange, "29", "Cranberry", "Strawberry")
TEMP_dataG29_StrawChe <- formattinglogchange(data_logchange, "29", "Strawberry", "Cherry")

2 ADAPTATION

2.1 Proportion of adaptation: raw data

tapply(data_sum$Nb_adults[data_sum$Generation=="2"],data_sum$Fruit_s[data_sum$Generation=="2"],mean)
## Blackcurrant       Cherry    Cranberry          Fig        Grape     Rosehips 
##           NA      8.80000     19.38667           NA           NA           NA 
##   Strawberry       Tomato 
##     12.19000           NA
tapply(data_sum$Nb_adults[data_sum$Generation=="27"],data_sum$Fruit_s[data_sum$Generation=="27"],mean)
## Blackcurrant       Cherry    Cranberry          Fig        Grape     Rosehips 
##           NA     35.84000     33.63750           NA           NA           NA 
##   Strawberry       Tomato 
##     42.33333           NA
## Calcul proportion change between G2 and G27 
((tapply(data_sum$Nb_adults[data_sum$Generation=="27"],data_sum$Fruit_s[data_sum$Generation=="27"],mean) -
    tapply(data_sum$Nb_adults[data_sum$Generation=="2"],data_sum$Fruit_s[data_sum$Generation=="2"],mean)) / 
  tapply(data_sum$Nb_adults[data_sum$Generation=="2"],data_sum$Fruit_s[data_sum$Generation=="2"],mean)) * 100
## Blackcurrant       Cherry    Cranberry          Fig        Grape     Rosehips 
##           NA    307.27273     73.50843           NA           NA           NA 
##   Strawberry       Tomato 
##    247.27919           NA
# 
# (abs(tapply(data_sum$fitness[data_sum$Generation=="27"],data_sum$Fruit_s[data_sum$Generation=="27"],mean) -
#     tapply(data_sum$fitness[data_sum$Generation=="2"],data_sum$Fruit_s[data_sum$Generation=="2"],mean)) / 
#   abs(tapply(data_sum$fitness[data_sum$Generation=="2"],data_sum$Fruit_s[data_sum$Generation=="2"],mean))) * 100


## Calcul proportion change between G2 and phase2 (because G8 is not representative)  
((tapply(data_sum$Nb_adults[data_sum$Generation=="8"],data_sum$Fruit_s[data_sum$Generation=="8"],mean) -
    tapply(data_sum$Nb_adults[data_sum$Generation=="2"],data_sum$Fruit_s[data_sum$Generation=="2"],mean)) / 
  tapply(data_sum$Nb_adults[data_sum$Generation=="2"],data_sum$Fruit_s[data_sum$Generation=="2"],mean)) * 100
## Blackcurrant       Cherry    Cranberry          Fig        Grape     Rosehips 
##           NA    289.61039     85.02183           NA           NA           NA 
##   Strawberry       Tomato 
##     89.19658           NA
## Compute proportion change between G2 and phase 2 (because G8 is not representative)  
((tapply(data_sum$Nb_adults[data_sum$Phase=="pool"],data_sum$Fruit_s[data_sum$Phase=="pool"],mean) -
    tapply(data_sum$Nb_adults[data_sum$Phase=="first_prepool"],data_sum$Fruit_s[data_sum$Phase=="first_prepool"],mean)) / 
  tapply(data_sum$Nb_adults[data_sum$Phase=="first_prepool"],data_sum$Fruit_s[data_sum$Phase=="first_prepool"],mean)) * 100
## Blackcurrant       Cherry    Cranberry          Fig        Grape     Rosehips 
##           NA    123.10237     93.62985           NA           NA           NA 
##   Strawberry       Tomato 
##     45.43089           NA
((tapply(data_sum$Nb_adults[data_sum$Phase=="second_postpool"],data_sum$Fruit_s[data_sum$Phase=="second_postpool"],mean) -
    tapply(data_sum$Nb_adults[data_sum$Phase=="first_prepool"],data_sum$Fruit_s[data_sum$Phase=="first_prepool"],mean)) / 
  tapply(data_sum$Nb_adults[data_sum$Phase=="first_prepool"],data_sum$Fruit_s[data_sum$Phase=="first_prepool"],mean)) * 100
## Blackcurrant       Cherry    Cranberry          Fig        Grape     Rosehips 
##           NA     230.1300     176.0760           NA           NA           NA 
##   Strawberry       Tomato 
##     145.1609           NA

2.2 Analysis

######## Models
mod1 <- lme4::lmer(fitness ~ 1 + (1|Generation) + (1|Generation:Fruit_s), 
            weights = N, data = data_sum, REML = FALSE) 
mod2 <- lme4::lmer(fitness~ Phase + (1|Generation) + (1|Generation:Fruit_s), 
            weights = N, data = data_sum, REML = FALSE)
mod3 <- lme4::lmer(fitness ~ Generation + (1|Generation) + (1|Generation:Fruit_s), 
            weights = N, data = data_sum, REML = FALSE) 
mod4 <- lme4::lmer(fitness ~ Fruit_s + (1|Generation) + (1|Generation:Fruit_s), 
            weights = N, data = data_sum, REML = FALSE) 
mod5 <- lme4::lmer(fitness ~ Fruit_s*Phase + (1|Generation) + (1|Generation:Fruit_s), 
            weights = N, data = data_sum, REML = FALSE) 
mod6 <- lme4::lmer(fitness ~ Fruit_s*Generation + (1|Generation) + (1|Generation:Fruit_s), 
            weights = N, data = data_sum, REML = FALSE) 
mod7 <- lme4::lmer(fitness ~ Fruit_s + Phase + (1|Generation) + (1|Generation:Fruit_s), 
            weights = N, data = data_sum, REML = FALSE) 
mod8 <- lme4::lmer(fitness ~ Fruit_s + Generation + (1|Generation) + (1|Generation:Fruit_s), 
            weights = N, data = data_sum, REML = FALSE) 

MuMIn::model.sel(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8)
## Model selection table 
##        (Int) Phs     Gnr Frt_s Frt_s:Phs Frt_s:Gnr df   logLik  AICc delta
## mod7 -0.3309   +             +                      8 -122.149 260.9  0.00
## mod2 -0.4770   +                                    6 -126.389 265.1  4.23
## mod5 -0.5293   +             +         +           12 -120.494 266.3  5.42
## mod8 -0.1564     0.03884     +                      7 -131.839 278.1 17.24
## mod6 -0.3216     0.04920     +                   +  9 -130.751 280.3 19.36
## mod3 -0.3093     0.03926                            5 -136.100 282.4 21.55
## mod4  0.4371                 +                      6 -139.518 291.4 30.48
## mod1  0.2875                                        4 -143.957 296.1 35.18
##      weight
## mod7  0.842
## mod2  0.102
## mod5  0.056
## mod8  0.000
## mod6  0.000
## mod3  0.000
## mod4  0.000
## mod1  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | Generation', '1 | Generation:Fruit_s'
summary(mod7)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: 
## fitness ~ Fruit_s + Phase + (1 | Generation) + (1 | Generation:Fruit_s)
##    Data: data_sum
## Weights: N
## 
##      AIC      BIC   logLik deviance df.resid 
##    260.3    288.4   -122.1    244.3      239 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2877 -0.3571  0.0899  0.4745  2.2266 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  Generation:Fruit_s (Intercept) 0.05857  0.2420  
##  Generation         (Intercept) 0.01257  0.1121  
##  Residual                       2.11342  1.4538  
## Number of obs: 247, groups:  Generation:Fruit_s, 72; Generation, 24
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)          -0.33091    0.11980  -2.762
## Fruit_sCranberry     -0.16282    0.08537  -1.907
## Fruit_sStrawberry    -0.26491    0.08767  -3.022
## Phasepool             0.57826    0.14805   3.906
## Phasesecond_postpool  0.99779    0.11899   8.386
## 
## Correlation of Fixed Effects:
##             (Intr) Frt_sC Frt_sS Phaspl
## Frt_sCrnbrr -0.385                     
## Frt_sStrwbr -0.381  0.532              
## Phasepool   -0.641 -0.025 -0.025       
## Phsscnd_pst -0.817  0.003  0.011  0.658
mod_Phase <- lme4::lmer(fitness ~ Fruit_s + Phase + (1|Generation) + (1|Generation:Fruit_s), 
            weights = N, data = data_sum) 


#Posthoc
emmeans::emmeans(mod_Phase, list(pairwise ~ Phase), adjust = "tukey") #
## $`emmeans of Phase`
##  Phase           emmean     SE    df lower.CL upper.CL
##  first_prepool   -0.474 0.1179   159   -0.758   -0.189
##  pool             0.106 0.1088 34492   -0.154    0.366
##  second_postpool  0.524 0.0549   543    0.392    0.655
## 
## Results are averaged over the levels of: Fruit_s 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Phase`
##  contrast                        estimate    SE   df t.ratio p.value
##  first_prepool - pool              -0.580 0.160  527 -3.617  0.0010 
##  first_prepool - second_postpool   -0.997 0.130  190 -7.671  <.0001 
##  pool - second_postpool            -0.418 0.122 6553 -3.429  0.0018 
## 
## Results are averaged over the levels of: Fruit_s 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates
emmeans::emmeans(mod_Phase, list(pairwise ~ Phase+Fruit_s), adjust = "tukey") #
## $`emmeans of Phase, Fruit_s`
##  Phase           Fruit_s     emmean     SE    df lower.CL upper.CL
##  first_prepool   Cherry     -0.3309 0.1296   232  -0.6926   0.0308
##  pool            Cherry      0.2490 0.1231 27016  -0.0915   0.5894
##  second_postpool Cherry      0.6665 0.0760  2263   0.4560   0.8770
##  first_prepool   Cranberry  -0.4941 0.1272   223  -0.8493  -0.1389
##  pool            Cranberry   0.0858 0.1179 32019  -0.2403   0.4118
##  second_postpool Cranberry   0.5033 0.0728   901   0.3015   0.7051
##  first_prepool   Strawberry -0.5961 0.1283   226  -0.9543  -0.2378
##  pool            Strawberry -0.0162 0.1190 34230  -0.3454   0.3130
##  second_postpool Strawberry  0.4014 0.0754  2177   0.1926   0.6101
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 9 estimates 
## 
## $`pairwise differences of Phase, Fruit_s`
##  contrast                                               estimate     SE    df
##  first_prepool,Cherry - pool,Cherry                       -0.580 0.1603   527
##  first_prepool,Cherry - second_postpool,Cherry            -0.997 0.1300   190
##  first_prepool,Cherry - first_prepool,Cranberry            0.163 0.0877  2089
##  first_prepool,Cherry - pool,Cranberry                    -0.417 0.1810   822
##  first_prepool,Cherry - second_postpool,Cranberry         -0.834 0.1573   365
##  first_prepool,Cherry - first_prepool,Strawberry           0.265 0.0898  3565
##  first_prepool,Cherry - pool,Strawberry                   -0.315 0.1820   869
##  first_prepool,Cherry - second_postpool,Strawberry        -0.732 0.1588   424
##  pool,Cherry - second_postpool,Cherry                     -0.418 0.1218  6553
##  pool,Cherry - first_prepool,Cranberry                     0.743 0.1845   936
##  pool,Cherry - pool,Cranberry                              0.163 0.0877  2089
##  pool,Cherry - second_postpool,Cranberry                  -0.254 0.1526  5336
##  pool,Cherry - first_prepool,Strawberry                    0.845 0.1855   953
##  pool,Cherry - pool,Strawberry                             0.265 0.0898  3565
##  pool,Cherry - second_postpool,Strawberry                 -0.152 0.1542  9464
##  second_postpool,Cherry - first_prepool,Cranberry          1.161 0.1565   449
##  second_postpool,Cherry - pool,Cranberry                   0.581 0.1475 13466
##  second_postpool,Cherry - second_postpool,Cranberry        0.163 0.0877  2089
##  second_postpool,Cherry - first_prepool,Strawberry         1.263 0.1573   422
##  second_postpool,Cherry - pool,Strawberry                  0.683 0.1484 10341
##  second_postpool,Cherry - second_postpool,Strawberry       0.265 0.0898  3565
##  first_prepool,Cranberry - pool,Cranberry                 -0.580 0.1603   527
##  first_prepool,Cranberry - second_postpool,Cranberry      -0.997 0.1300   190
##  first_prepool,Cranberry - first_prepool,Strawberry        0.102 0.0861  2174
##  first_prepool,Cranberry - pool,Strawberry                -0.478 0.1819   882
##  first_prepool,Cranberry - second_postpool,Strawberry     -0.895 0.1563   440
##  pool,Cranberry - second_postpool,Cranberry               -0.418 0.1218  6553
##  pool,Cranberry - first_prepool,Strawberry                 0.682 0.1820   848
##  pool,Cranberry - pool,Strawberry                          0.102 0.0861  2174
##  pool,Cranberry - second_postpool,Strawberry              -0.316 0.1495 13058
##  second_postpool,Cranberry - first_prepool,Strawberry      1.099 0.1556   356
##  second_postpool,Cranberry - pool,Strawberry               0.520 0.1487  5535
##  second_postpool,Cranberry - second_postpool,Strawberry    0.102 0.0861  2174
##  first_prepool,Strawberry - pool,Strawberry               -0.580 0.1603   527
##  first_prepool,Strawberry - second_postpool,Strawberry    -0.997 0.1300   190
##  pool,Strawberry - second_postpool,Strawberry             -0.418 0.1218  6553
##  t.ratio p.value
##  -3.617  0.0098 
##  -7.671  <.0001 
##   1.860  0.6416 
##  -2.302  0.3420 
##  -5.304  <.0001 
##   2.952  0.0772 
##  -1.729  0.7284 
##  -4.612  0.0002 
##  -3.429  0.0177 
##   4.027  0.0020 
##   1.860  0.6416 
##  -1.666  0.7672 
##   4.554  0.0002 
##   2.952  0.0772 
##  -0.988  0.9870 
##   7.418  <.0001 
##   3.937  0.0027 
##   1.860  0.6416 
##   8.027  <.0001 
##   4.602  0.0001 
##   2.952  0.0772 
##  -3.617  0.0098 
##  -7.671  <.0001 
##   1.184  0.9599 
##  -2.627  0.1770 
##  -5.729  <.0001 
##  -3.429  0.0177 
##   3.746  0.0060 
##   1.184  0.9599 
##  -2.110  0.4661 
##   7.065  <.0001 
##   3.493  0.0142 
##   1.184  0.9599 
##  -3.617  0.0098 
##  -7.671  <.0001 
##  -3.429  0.0177 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 9 estimates

2.3 Plot

#Position dodge
pd <- ggplot2::position_dodge(0.3) # move them .05 to the left and right

#Extract slope and intercept
dat_predict_allfruits <- expand.grid(Generation=as.numeric(levels(as.factor(data_sum$Generation))),
                                     Fruit_s=unique(data_sum$Fruit_s))

dat_predict_allfruits$Phase <- ifelse(dat_predict_allfruits$Generation<= 7, "first_prepool",
                                     ifelse(dat_predict_allfruits$Generation>=12, "second_postpool", "pool"))
dat_predict_allfruits$fitness_predicted <- predict(mod_Phase, newdata = dat_predict_allfruits, 
                          re.form= NA, type = "response")

#REAL DATA
#Add G1 G6 and G7
TEMP_lineG1 <- c(rep(NA, 2), 1,rep(NA, 5))
TEMP_lineG6 <- c(rep(NA, 2), 6,rep(NA, 5))
TEMP_lineG7 <- c(rep(NA, 2), 7,rep(NA, 5))

TEMP_total <- rbind(data_sum,TEMP_lineG1,TEMP_lineG1,TEMP_lineG1,TEMP_lineG6,TEMP_lineG6,TEMP_lineG6,TEMP_lineG7,TEMP_lineG7,TEMP_lineG7)
TEMP_total$Fruit_s[TEMP_total$Generation == "6"] <- c("Strawberry", "Cranberry", "Cherry")
TEMP_total$Fruit_s[TEMP_total$Generation == "1"] <- c("Strawberry", "Cranberry", "Cherry")
TEMP_total$Fruit_s[TEMP_total$Generation == "7"] <- c("Strawberry", "Cranberry", "Cherry")
tail(TEMP_total)
##     Line    Fruit_s Generation Phase  N Nb_adults sd fitness se_fitness
## 251 <NA> Strawberry          6  <NA> NA        NA NA      NA         NA
## 252 <NA>  Cranberry          6  <NA> NA        NA NA      NA         NA
## 253 <NA>     Cherry          6  <NA> NA        NA NA      NA         NA
## 254 <NA> Strawberry          7  <NA> NA        NA NA      NA         NA
## 255 <NA>  Cranberry          7  <NA> NA        NA NA      NA         NA
## 256 <NA>     Cherry          7  <NA> NA        NA NA      NA         NA
## Add label 
TEMP_anno <- data.frame(x1 = c(3.5, 3.5,  10, 3.5,  3.5, 10,  3.5, 3.5, 10), 
                        x2 = c(9, 17.5, 17.5,9, 17.5,  17.5,9,17.5, 17.5),
                        y1 = c(1,2,1.2,1,2,1.2,1,2,1.2), 
                        y2 = c(1.25,  2.25, 1.45, 1.25,  2.25, 1.45,1.25,  2.25, 1.45),
                        xstar = c(6.5,10,14,6.5,10,14,6.5,10,14), 
                        ystar = c(1.5,2.5,1.7,1.5,2.5,1.7,1.5,2.5,1.7),
                        lab = c("**", "***", "**", "**",  "***",  "**", "**" ,"***",  "**"),
                        Fruit_s = c("Cherry", "Cherry", "Cherry",
                               "Cranberry", "Cranberry", "Cranberry",
                               "Strawberry","Strawberry","Strawberry"), 
                        Line = NA)
TEMP_anno
##     x1   x2  y1   y2 xstar ystar lab    Fruit_s Line
## 1  3.5  9.0 1.0 1.25   6.5   1.5  **     Cherry   NA
## 2  3.5 17.5 2.0 2.25  10.0   2.5 ***     Cherry   NA
## 3 10.0 17.5 1.2 1.45  14.0   1.7  **     Cherry   NA
## 4  3.5  9.0 1.0 1.25   6.5   1.5  **  Cranberry   NA
## 5  3.5 17.5 2.0 2.25  10.0   2.5 ***  Cranberry   NA
## 6 10.0 17.5 1.2 1.45  14.0   1.7  **  Cranberry   NA
## 7  3.5  9.0 1.0 1.25   6.5   1.5  ** Strawberry   NA
## 8  3.5 17.5 2.0 2.25  10.0   2.5 *** Strawberry   NA
## 9 10.0 17.5 1.2 1.45  14.0   1.7  ** Strawberry   NA
TEMP_title <- data.frame(xtitle = c(14, 14, 14),
                       ytitle = c(3, 3, 3),
                       title = c("Cherry", "Cranberry", "Strawberry"),
                       Fruit_s = c("Cherry", "Cranberry", "Strawberry"), 
                       Line = NA)
TEMP_title
##   xtitle ytitle      title    Fruit_s Line
## 1     14      3     Cherry     Cherry   NA
## 2     14      3  Cranberry  Cranberry   NA
## 3     14      3 Strawberry Strawberry   NA
PLOT_FITNESS_CHERRY <- ggplot(data = TEMP_total[TEMP_total$Fruit_s == "Cherry",], 
                            aes(x = factor(Generation),group = Line, y = fitness, colour =Fruit_s)) + 
  geom_errorbar(aes(ymin =fitness-1.96*se_fitness, ymax = fitness + 1.96*se_fitness),
                width=.1,position = pd, size = 0.2,color = "black") + 
  geom_line(size = 0.3,position = pd) + 
  geom_line(data = dat_predict_allfruits[dat_predict_allfruits$Fruit_s == "Cherry",], 
                     aes(x = factor(Generation), y = fitness_predicted,
                         colour = "black", group = Phase), size = 0.5) + 
  geom_point(size =1, position = pd, shape =21, fill = "white") + 
  ylim(-3, 3.05) + 
  ylab("Fitness") + 
  xlab("Generation") + 
  geom_text(data = TEMP_anno[TEMP_anno$Fruit_s == "Cherry",], aes(x = xstar,  y = ystar, label = lab), size =3.3) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cherry",], aes(x = x1, xend = x1, 
           y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cherry",], aes(x = x2, xend = x2, 
           y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cherry",], aes(x = x1, xend = x2, 
           y = y2, yend = y2), size = 0.4) + 
  scale_color_manual(values = c("black", "#BC3C6D", "#FDB424", "#3FAA96")) + 
  ggtitle("Cherry") + 
  theme_LO_adaptation + theme(plot.title = element_text(color = "#BC3C6D"))
PLOT_FITNESS_CHERRY

PLOT_FITNESS_CRANB <- ggplot2::ggplot(data = TEMP_total[TEMP_total$Fruit_s == "Cranberry",], 
                            aes(x = factor(Generation),group = Line, y = fitness, colour =Fruit_s)) + 
  geom_errorbar(aes(ymin =fitness-1.96*se_fitness, ymax = fitness + 1.96*se_fitness),
                width=.1,position = pd, size = 0.2,color = "black") + 
  geom_line(size = 0.3,position = pd) + 
  geom_line(data = dat_predict_allfruits[dat_predict_allfruits$Fruit_s == "Cranberry",], aes(x = factor(Generation), y = fitness_predicted, colour = "black", group = Phase),
                size = 0.5) + 
  geom_point(size =1, position = pd, shape =21, fill = "white") + 
  ylim(-3, 3.05) + 
  ylab("Fitness") + 
  xlab("Generation") + 
  geom_text(data = TEMP_anno[TEMP_anno$Fruit_s == "Cranberry",], aes(x = xstar,  y = ystar, label = lab), size =3.3) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cranberry",], aes(x = x1, xend = x1, 
           y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cranberry",], aes(x = x2, xend = x2, 
           y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Cranberry",], aes(x = x1, xend = x2, 
           y = y2, yend = y2), size = 0.4) + 
  scale_color_manual(values = c("black", "#FDB424")) + 
  ggtitle("Cranberry") + 
  theme_LO_adaptation + theme(plot.title = element_text(color = "#FDB424"))
PLOT_FITNESS_CRANB

PLOT_FITNESS_STRAW <- ggplot2::ggplot(data = TEMP_total[TEMP_total$Fruit_s == "Strawberry",], 
                            aes(x = factor(Generation),group = Line, y = fitness, colour =Fruit_s)) + 
  geom_errorbar(aes(ymin =fitness-1.96*se_fitness, ymax = fitness + 1.96*se_fitness),
                width=.1,position = pd, size = 0.2,color = "black") + 
  geom_line(size = 0.3,position = pd) + 
  geom_line(data = dat_predict_allfruits[dat_predict_allfruits$Fruit_s == "Strawberry",], 
            aes(x = factor(Generation), y = fitness_predicted, colour = "black", group = Phase),
            size = 0.5) + 
  geom_point(size =1, position = pd, shape =21, fill = "white") + 
  ylim(-3, 3.05) + 
  ylab("Fitness") + 
  xlab("Generation") + 
  geom_text(data = TEMP_anno[TEMP_anno$Fruit_s == "Strawberry",], 
            aes(x = xstar,  y = ystar, label = lab), size =3.3) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Strawberry",], 
               aes(x = x1, xend = x1, y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Strawberry",], 
               aes(x = x2, xend = x2, y = y1, yend = y2), size = 0.4) + 
  geom_segment(data = TEMP_anno[TEMP_anno$Fruit_s == "Strawberry",], 
               aes(x = x1, xend = x2, y = y2, yend = y2), size = 0.4) + 
  scale_color_manual(values = c("black", "#3FAA96")) + 
  ggtitle("Strawberry") + 
  theme_LO_adaptation + theme(plot.title = element_text(color = "#3FAA96"))
PLOT_FITNESS_STRAW

DYNAMIQUE_THREE_JOIN <- cowplot::ggdraw() + 
  cowplot::draw_plot(PLOT_FITNESS_CHERRY + theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank()), 
            x = 0, y = 0.66, width = 1, height = 0.33) + 
  cowplot::draw_plot(PLOT_FITNESS_CRANB + theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.line.x = element_blank()), 
            x = 0, y = 0.33, width = 1, height = 0.33) + 
  cowplot::draw_plot(PLOT_FITNESS_STRAW, 
            x = 0, y = 00, width = 1, height = 0.33) + 
  cowplot::draw_plot_label(c("A", "B", "C"),  
                  x = c(0, 0, 0), 
                  y = c(1, 0.66, 0.33), 
                  hjust = c(-0.5, -0.5, -0.5), 
                  vjust = c(1.5, 1.5, 1.5),
                  size = 12) 
 DYNAMIQUE_THREE_JOIN

cowplot::save_plot(file =here::here("figures", "FIG_Adaptation.pdf"), DYNAMIQUE_THREE_JOIN, base_height = 17/cm(1), base_width = 11/cm(1), dpi = 1200)

3 HETEROGENEITY

3.1 Plot fitness

pd <-  position_dodge(width = 0.5)



################### INTERMEDIATE PHENOTYPING

# Re-order levels of Line
data_sum_G7 <- data_logchange[data_logchange$Generation == "7",]
data_sum_G7$Line <- factor(data_sum_G7$Line, levels= c("CE2", "CE1", "CE4", "CE3",
                                                       "CR2", "CR3", "CR5", "CR1", "CR4",
                                                       "FR2", "FR3", "FR5", "FR1", "FR4"))


#Plot: 
symp_allop_g7 <- ggplot(data = data_sum_G7, 
                      aes(x = Line, y = logchange, group = Treatment, 
                          color = Fruit_s, shape = Treatment, fill = Fruit_s)) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size = 0.5) + 
  geom_errorbar(aes(ymin = logchange-1.96*sd_logchange, 
                    ymax =logchange + 1.96*sd_logchange, linetype = Line_type),
                  width = 0.2, size = 0.5, alpha = 0.6,position = pd) + 
  geom_point(position = pd, fill = "white", size =3) + 
  ylab("Fitness change between intermediate\nand initial phenotyping steps")  + 
  xlab("Populations") + 
  labs(shape = "Test fruit", color = "Selection fruit") + 
  guides(color = FALSE,
           shape = guide_legend(override.aes = list(fill = c("black")))) + 
  scale_shape_manual(values = c(21, 22, 24)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
  scale_fill_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
  theme_LO_sober + theme(axis.text.x  = element_blank()) +
  coord_cartesian(expand = FALSE, ylim = c(-2, 2), xlim=c(0, 15),clip = "off") + 
  ggtitle("Intermediate phenotyping step") + 
  annotate("text", x = 2.5, y =1.5, label = 'bold(" Evolved\non cherry")',
                     size =3,colour = "#BC3C6D",parse = TRUE) + 
  annotate("text", x = 7.5, y =1.5, label = 'bold("   Evolved\non cranberry")',
                     size =3,colour = "#FDB424",parse = TRUE) + 
  annotate("text", x = 12.5, y =1.5, label = 'bold("   Evolved\non strawberry")',
                     size =3,colour = "#3FAA96",parse = TRUE) + 
  geom_segment(x = 15.5, y = -0.1, xend= 15.5, yend = -1.2, size = 0.15, 
               arrow = arrow(length = unit(0.05, "npc")),colour = "black") + 
  geom_segment(x = 15.5, y = 0.1, xend= 15.5, yend = 1.2, size = 0.15,
               arrow = arrow(length = unit(0.05, "npc")),colour = "black") + 
  annotate("text", x = 16.2, y = 0.5, label = 'bold(" Fitness\nincrease")',
             size =3,colour = "black",parse = TRUE, angle =90) + 
  annotate("text", x = 16.2, y =-0.5, label = 'bold("  Fitness\ndecrease")',
              size =3,colour = "black",parse = TRUE, angle =90) 
symp_allop_g7

################### FINAL PHENOTYPING

# Re-order levels of Line
data_sum_G29 <- data_logchange[data_logchange$Generation == "29",]
data_sum_G29$Line <- factor(data_sum_G29$Line, levels = c("CEA", "CEB", "CEC", 
                                                          "CRD", "CRA", "CRC", "CRB", "CRE",
                                                          "FRA", "FRC", "FRB"))

symp_allop_G29 <- ggplot(data = data_sum_G29, 
              aes(x = Line, y=logchange, group = Treatment, 
                  color =Fruit_s, shape = Treatment,fill =Fruit_s)) + 
      geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size = 0.5) + 
    geom_point(size =3, position = pd) + 
    geom_errorbar(aes(ymin =logchange-1.96*sd_logchange, ymax =logchange + 1.96*sd_logchange),
                width= 0.2, size = 0.5, alpha = 0.6,position = pd) + 
    ylab("Fitness change between final\nand initial phenotyping steps")  + 
    xlab("Populations") + 
    labs(shape = "Test fruit", color = "Selection fruit") + 
    guides(fill = FALSE, alpha = FALSE) + 
    scale_shape_manual(values = c(21, 22, 24)) + 
    scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
    scale_fill_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
    theme_LO_sober + theme(axis.text.x  = element_blank()) +
    coord_cartesian(expand = FALSE, ylim = c(-2, 2), xlim=c(0, 11),clip = "off") + 
  annotate("text", x = 2, y =1.5, label = 'bold(" Evolved\non cherry")',
                     size =3,colour = "#BC3C6D",parse = TRUE) + 
  annotate("text", x = 6, y =1.5, label = 'bold("   Evolved\non cranberry")',
                     size =3,colour = "#FDB424",parse = TRUE) + 
  annotate("text", x = 10, y =1.5, label = 'bold("   Evolved\non strawberry")',
                     size =3,colour = "#3FAA96",parse = TRUE) + 
  geom_segment(x = 11.5, y = -0.1, xend= 11.5, yend = -1.2, size = 0.15, 
               arrow = arrow(length = unit(0.05, "npc")),colour = "black") + 
    geom_segment(x = 11.5, y = 0.1, xend= 11.5, yend = 1.2, size = 0.15,
               arrow = arrow(length = unit(0.05, "npc")),colour = "black") + 
    annotate("text", x = 12, y = 0.5, label = 'bold(" Fitness\nincrease")',
             size =3,colour = "black",parse = TRUE, angle =90) + 
    annotate("text", x = 12, y = -0.5, label = 'bold("  Fitness\ndecrease")',
              size =3,colour = "black", parse = TRUE, angle =90) + 
  ggtitle("Final phenotyping step") 
symp_allop_G29

### Add stroke
DIF_FITNESS_G7 <- symp_allop_g7 + 
  geom_point(aes(alpha = SA, color = interaction(SA,Fruit_s)), position = pd, size =3, stroke =1.5, fill = "white") + 
  scale_alpha_manual(values = c(0, 1)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96", "#7C2748", "#CA8702", "#328677", "#BC3C6D", "#FDB424", "#3FAA96"))  
DIF_FITNESS_G7

DIF_FITNESS_G29 <- symp_allop_G29 + geom_point(aes(alpha = SA, color = interaction(SA,Fruit_s)), position = pd, size =3, stroke =1.5) + 
      scale_alpha_manual(values = c(0, 1)) + 
      scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96", "#7C2748", "#CA8702", "#328677", "#BC3C6D", "#FDB424", "#3FAA96"))  
DIF_FITNESS_G29

plot_legend <- ggplot(data = data_logchange, aes(x = Line, y=logchange, group = Treatment, 
                  color =Fruit_s, shape = Treatment,fill =Fruit_s)) + 
    geom_point(size =3, position = pd) + 
    labs(shape = "Test fruit", color = "Selection fruit") + 
    scale_shape_manual(values = c(21, 22, 24)) + 
    scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
    theme_LO_sober + 
    guides(fill = FALSE) + 
    theme(legend.key.size = unit(0.5, "cm"),
        axis.text.x  = element_blank()) +
    coord_cartesian(expand = FALSE, ylim = c(-2, 2), xlim=c(0, 11),clip = "off") 


###########################################################
#############################################ALL PLOT
legend_trade <- lemon::g_legend(plot_legend)
    

SYMP_ALLOP_TOTAL <- cowplot::ggdraw() + 
  cowplot::draw_plot(DIF_FITNESS_G7 + theme(legend.position = 'none', 
                             plot.margin =unit(c(0.5, 3.5, 0.5, 0.5), "cm")),
            x = 0, y = 0.5, width = 1, height = 0.5) + 
  cowplot::draw_plot(DIF_FITNESS_G29 + theme(legend.position = 'none', 
                             plot.margin =unit(c(0.5, 3.5, 0.5, 0.5), "cm")), 
            x = 0, y = 0, width = 1, height = 0.5) + 
  cowplot::draw_plot_label(c("A", "B"),  
                  x = c(0, 0), y = c(1, 0.5), 
                  size = 14) + 
  cowplot::draw_plot(legend_trade, x = 0.93, y = 0.5, width = 0.001, height = 0.001) 

SYMP_ALLOP_TOTAL

cowplot::save_plot(file =here::here("figures", "FIG_Heterogeneity.pdf"),  
                   SYMP_ALLOP_TOTAL, base_height = 20/cm(1), base_width = 20/cm(1), dpi = 1200)

4 Traits

4.1 Correlation between Nb_eggs and emergence_rate

tapply(data_G0$Nb_eggs, data_G0$Treatment, mean)
##     Cherry  Cranberry Strawberry 
##     149.49     122.85     107.84
diag(tapply(data_G7$Nb_eggs, list(data_G7$Treatment,data_G7$Fruit_s), mean))
##     Cherry  Cranberry Strawberry 
##   195.4375   181.1591   174.2923
diag(tapply(data_G29$Nb_eggs, list(data_G29$Treatment,data_G29$Fruit_s), mean))
##     Cherry  Cranberry Strawberry 
##   103.7667   123.7000   126.7333
tapply(data_G0$Emergence_rate, data_G0$Treatment, mean)
##     Cherry  Cranberry Strawberry 
##  0.1837216  0.2660867  0.2777895
diag(tapply(data_G7$Emergence_rate, list(data_G7$Treatment,data_G7$Fruit_s), mean))
##     Cherry  Cranberry Strawberry 
##  0.1816565  0.2070999  0.2728989
diag(tapply(data_G29$Emergence_rate, list(data_G29$Treatment,data_G29$Fruit_s), mean))
##     Cherry  Cranberry Strawberry 
##  0.3421827  0.2794123  0.3543303
#Plot

Nb_eggsEmergence_rateG7 <- ggplot(data = data_logchangeNb_eggsEmergence_rate[data_logchangeNb_eggsEmergence_rate$Generation==7,],
                          aes(x = Nb_eggslogchange,y = Emergence_ratelogchange, group = Treatment, shape = Treatment)) +
  geom_errorbarh(aes(xmin = Nb_eggslogchange-1.96*sdNb_eggslogchange,
                     xmax = Nb_eggslogchange+1.96*sdNb_eggslogchange),
                 height=0.02,size=0.2,alpha=1) +
  geom_errorbar(aes(ymin = Emergence_ratelogchange-1.96*sdEmergence_ratelogchange,
                    ymax = Emergence_ratelogchange+1.96*sdEmergence_ratelogchange),
                width=0.02,size=0.2,alpha=1) + 
  geom_point(aes(shape=Treatment, color=Treatment)) +
  
  xlab("Change in number of eggs between\nintermediate and initial phenotyping steps")  + 
  ylab("Change in egg-to-adult between\nintermediate and initial phenotyping steps")  +
  xlim(-0.75, 1) + 
  ylim(-1.75, 1.5) +   
  labs(shape = "Test fruit") + 
  scale_shape_manual(values = c(21, 22, 24)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
  scale_fill_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
  theme_LO_sober
 Nb_eggsEmergence_rateG7 <- Nb_eggsEmergence_rateG7 + facet_wrap(~ Fruit_s)

 
 
 Nb_eggsEmergence_rateG29 <- ggplot(data = data_logchangeNb_eggsEmergence_rate[data_logchangeNb_eggsEmergence_rate$Generation==29,],
                          aes(x = Nb_eggslogchange,y = Emergence_ratelogchange, group = Treatment, shape = Treatment)) +
  geom_errorbarh(aes(xmin = Nb_eggslogchange-1.96*sdNb_eggslogchange,
                     xmax = Nb_eggslogchange+1.96*sdNb_eggslogchange),
                 height=0.02,size=0.2,alpha=1) +
  geom_errorbar(aes(ymin = Emergence_ratelogchange-1.96*sdEmergence_ratelogchange,
                    ymax = Emergence_ratelogchange+1.96*sdEmergence_ratelogchange),
                width=0.02,size=0.2,alpha=1) + 
  geom_point(aes(shape=Treatment, color=Treatment)) +
  
  xlab("Change in number of eggs between\nfinal and initial phenotyping steps")  + 
  ylab("Change in egg-to-adult viability between\nfinal and initial phenotyping steps")  +
  xlim(-0.75, 1) + 
  ylim(-1.75, 1.5) +   
  labs(shape = "Test fruit", color = "Test fruit") + 
  scale_shape_manual(values = c(21, 22, 24)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
  scale_fill_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
  theme_LO_sober
Nb_eggsEmergence_rateG29 <- Nb_eggsEmergence_rateG29 + facet_wrap(~ Fruit_s)
 
legend_trade <- lemon::g_legend(Nb_eggsEmergence_rateG29)

Nb_eggsEmergence_rate <- cowplot::ggdraw() + 
  cowplot::draw_plot(Nb_eggsEmergence_rateG7 + theme(legend.position = 'none', 
                             plot.margin =unit(c(0.5, 3.5, 0.5, 0.5), "cm")),
            x = 0, y = 0.5, width = 1, height = 0.5) + 
  cowplot::draw_plot(Nb_eggsEmergence_rateG29 + theme(legend.position = 'none', 
                             plot.margin =unit(c(0.5, 3.5, 0.5, 0.5), "cm")), 
            x = 0, y = 0, width = 1, height = 0.5) + 
  cowplot::draw_plot_label(c("A", "B"),  
                  x = c(0, 0), y = c(1, 0.5), 
                  size = 14) + 
  cowplot::draw_plot(legend_trade, x = 0.93, y = 0.5, width = 0.001, height = 0.001) 

cowplot::save_plot(file =here::here("figures", "FIG_SX_Nb_eggsEmergence_rate.pdf"),  
                   Nb_eggsEmergence_rate, base_height = 20/cm(1), base_width = 20/cm(1), dpi = 1200)

4.2 Plot Nb_eggs

pd <-  position_dodge(width = 0.5)



################### INTERMEDIATE PHENOTYPING

# Re-order levels of Line
data_sum_G7 <- data_logchangeNb_eggs[data_logchangeNb_eggs$Generation == "7",]
data_sum_G7$Line <- factor(data_sum_G7$Line, levels= c("CE2", "CE1", "CE4", "CE3",
                                                       "CR2", "CR3", "CR5", "CR1", "CR4",
                                                       "FR2", "FR3", "FR5", "FR1", "FR4"))

#Plot: 
symp_allop_g7 <- ggplot(data = data_sum_G7, 
                      aes(x = Line, y = logchange, group = Treatment, 
                          color = Fruit_s, shape = Treatment, fill = Fruit_s)) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size = 0.5) + 
  geom_errorbar(aes(ymin = logchange-1.96*sd_logchange, 
                    ymax =logchange + 1.96*sd_logchange, linetype = Line_type),
                  width = 0.2, size = 0.5, alpha = 0.6,position = pd) + 
  geom_point(position = pd, fill = "white", size =3) + 
  ylab("Change in number of eggs between\nintermediate and initial phenotyping steps")  + 
  xlab("Populations") + 
  labs(shape = "Test fruit", color = "Selection fruit") + 
  guides(color = FALSE,
           shape = guide_legend(override.aes = list(fill = c("black")))) + 
  scale_shape_manual(values = c(21, 22, 24)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
  scale_fill_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
  theme_LO_sober + theme(axis.text.x  = element_blank()) +
  coord_cartesian(expand = FALSE, ylim = c(-2, 2), xlim=c(0, 15),clip = "off") + 
  ggtitle("Intermediate phenotyping step") + 
  annotate("text", x = 2.5, y =1.5, label = 'bold(" Evolved\non cherry")',
                     size =3,colour = "#BC3C6D",parse = TRUE) + 
  annotate("text", x = 7.5, y =1.5, label = 'bold("   Evolved\non cranberry")',
                     size =3,colour = "#FDB424",parse = TRUE) + 
  annotate("text", x = 12.5, y =1.5, label = 'bold("   Evolved\non strawberry")',
                     size =3,colour = "#3FAA96",parse = TRUE)
symp_allop_g7

################### FINAL PHENOTYPING

# Re-order levels of Line
data_sum_G29 <- data_logchangeNb_eggs[data_logchangeNb_eggs$Generation == "29",]
data_sum_G29$Line <- factor(data_sum_G29$Line, levels = c("CEA", "CEB", "CEC", 
                                                          "CRD", "CRA", "CRC", "CRB", "CRE",
                                                          "FRA", "FRC", "FRB"))

symp_allop_G29 <- ggplot(data = data_sum_G29, 
              aes(x = Line, y=logchange, group = Treatment, 
                  color =Fruit_s, shape = Treatment,fill =Fruit_s)) + 
      geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size = 0.5) + 
    geom_point(size =3, position = pd) + 
    geom_errorbar(aes(ymin =logchange-1.96*sd_logchange, ymax =logchange + 1.96*sd_logchange),
                width= 0.2, size = 0.5, alpha = 0.6,position = pd) + 
    ylab("Change in number of eggs between\nfinal and initial phenotyping steps")  + 
    xlab("Populations") + 
    labs(shape = "Test fruit", color = "Selection fruit") + 
    guides(fill = FALSE, alpha = FALSE) + 
    scale_shape_manual(values = c(21, 22, 24)) + 
    scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
    scale_fill_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
    theme_LO_sober + theme(axis.text.x  = element_blank()) +
    coord_cartesian(expand = FALSE, ylim = c(-2, 2), xlim=c(0, 11),clip = "off") + 
  annotate("text", x = 2, y =1.5, label = 'bold(" Evolved\non cherry")',
                     size =3,colour = "#BC3C6D",parse = TRUE) + 
  annotate("text", x = 6, y =1.5, label = 'bold("   Evolved\non cranberry")',
                     size =3,colour = "#FDB424",parse = TRUE) + 
  annotate("text", x = 10, y =1.5, label = 'bold("   Evolved\non strawberry")',
                     size =3,colour = "#3FAA96",parse = TRUE) + 
  ggtitle("Final phenotyping step") 
symp_allop_G29

### Add stroke
DIF_FITNESS_G7 <- symp_allop_g7 + 
  geom_point(aes(alpha = SA, color = interaction(SA,Fruit_s)), position = pd, size =3, stroke =1.5, fill = "white") + 
  scale_alpha_manual(values = c(0, 1)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96", "#7C2748", "#CA8702", "#328677", "#BC3C6D", "#FDB424", "#3FAA96"))  
DIF_FITNESS_G7

DIF_FITNESS_G29 <- symp_allop_G29 + geom_point(aes(alpha = SA, color = interaction(SA,Fruit_s)), position = pd, size =3, stroke =1.5) + 
      scale_alpha_manual(values = c(0, 1)) + 
      scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96", "#7C2748", "#CA8702", "#328677", "#BC3C6D", "#FDB424", "#3FAA96"))  
DIF_FITNESS_G29

plot_legend <- ggplot(data = data_logchange, aes(x = Line, y=logchange, group = Treatment, 
                  color =Fruit_s, shape = Treatment,fill =Fruit_s)) + 
    geom_point(size =3, position = pd) + 
    labs(shape = "Test fruit", color = "Selection fruit") + 
    scale_shape_manual(values = c(21, 22, 24)) + 
    scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
    theme_LO_sober + 
    guides(fill = FALSE) + 
    theme(legend.key.size = unit(0.5, "cm"),
        axis.text.x  = element_blank()) +
    coord_cartesian(expand = FALSE, ylim = c(-2, 2), xlim=c(0, 11),clip = "off") 


###########################################################
#############################################ALL PLOT
legend_trade <- lemon::g_legend(plot_legend)
    

SYMP_ALLOP_TOTAL <- cowplot::ggdraw() + 
  cowplot::draw_plot(DIF_FITNESS_G7 + theme(legend.position = 'none', 
                             plot.margin =unit(c(0.5, 3.5, 0.5, 0.5), "cm")),
            x = 0, y = 0.5, width = 1, height = 0.5) + 
  cowplot::draw_plot(DIF_FITNESS_G29 + theme(legend.position = 'none', 
                             plot.margin =unit(c(0.5, 3.5, 0.5, 0.5), "cm")), 
            x = 0, y = 0, width = 1, height = 0.5) + 
  cowplot::draw_plot_label(c("A", "B"),  
                  x = c(0, 0), y = c(1, 0.5), 
                  size = 14) + 
  cowplot::draw_plot(legend_trade, x = 0.93, y = 0.5, width = 0.001, height = 0.001) 

DIF_FITNESS_G7Nb_eggs <- DIF_FITNESS_G7
DIF_FITNESS_G29Nb_eggs <- DIF_FITNESS_G29

cowplot::save_plot(file =here::here("figures", "FIG_SX_HeterogeneityNb_eggs.pdf"),  
                   SYMP_ALLOP_TOTAL, base_height = 20/cm(1), base_width = 20/cm(1), dpi = 1200)

4.3 Plot egg-to-adult viability

pd <-  position_dodge(width = 0.5)



################### INTERMEDIATE PHENOTYPING

# Re-order levels of Line
data_sum_G7 <- data_logchangeEmergence_rate[data_logchangeEmergence_rate$Generation == "7",]
data_sum_G7$Line <- factor(data_sum_G7$Line, levels= c("CE2", "CE1", "CE4", "CE3",
                                                       "CR2", "CR3", "CR5", "CR1", "CR4",
                                                       "FR2", "FR3", "FR5", "FR1", "FR4"))

#Plot: 
symp_allop_g7 <- ggplot(data = data_sum_G7, 
                      aes(x = Line, y = logchange, group = Treatment, 
                          color = Fruit_s, shape = Treatment, fill = Fruit_s)) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size = 0.5) + 
  geom_errorbar(aes(ymin = logchange-1.96*sd_logchange, 
                    ymax =logchange + 1.96*sd_logchange, linetype = Line_type),
                  width = 0.2, size = 0.5, alpha = 0.6,position = pd) + 
  geom_point(position = pd, fill = "white", size =3) + 
  ylab("Change in egg-to-adult viability\nbetween intermediate and initial phenotyping steps")  + 
  xlab("Populations") + 
  labs(shape = "Test fruit", color = "Selection fruit") + 
  guides(color = FALSE,
           shape = guide_legend(override.aes = list(fill = c("black")))) + 
  scale_shape_manual(values = c(21, 22, 24)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
  scale_fill_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
  theme_LO_sober + theme(axis.text.x  = element_blank()) +
  coord_cartesian(expand = FALSE, ylim = c(-2, 2), xlim=c(0, 15),clip = "off") + 
  ggtitle("Intermediate phenotyping step") + 
  annotate("text", x = 2.5, y =1.5, label = 'bold(" Evolved\non cherry")',
                     size =3,colour = "#BC3C6D",parse = TRUE) + 
  annotate("text", x = 7.5, y =1.5, label = 'bold("   Evolved\non cranberry")',
                     size =3,colour = "#FDB424",parse = TRUE) + 
  annotate("text", x = 12.5, y =1.5, label = 'bold("   Evolved\non strawberry")',
                     size =3,colour = "#3FAA96",parse = TRUE)
symp_allop_g7

################### FINAL PHENOTYPING

# Re-order levels of Line
data_sum_G29 <- data_logchangeEmergence_rate[data_logchangeEmergence_rate$Generation == "29",]
data_sum_G29$Line <- factor(data_sum_G29$Line, levels = c("CEA", "CEB", "CEC", 
                                                          "CRD", "CRA", "CRC", "CRB", "CRE",
                                                          "FRA", "FRC", "FRB"))

symp_allop_G29 <- ggplot(data = data_sum_G29, 
              aes(x = Line, y=logchange, group = Treatment, 
                  color =Fruit_s, shape = Treatment,fill =Fruit_s)) + 
      geom_hline(yintercept = 0, linetype = "dashed", color = "grey", size = 0.5) + 
    geom_point(size =3, position = pd) + 
    geom_errorbar(aes(ymin =logchange-1.96*sd_logchange, ymax =logchange + 1.96*sd_logchange),
                width= 0.2, size = 0.5, alpha = 0.6,position = pd) + 
    ylab("Change in egg-to-adult viability\nbetween final and initial phenotyping steps")  + 
    xlab("Populations") + 
    labs(shape = "Test fruit", color = "Selection fruit") + 
    guides(fill = FALSE, alpha = FALSE) + 
    scale_shape_manual(values = c(21, 22, 24)) + 
    scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
    scale_fill_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
    theme_LO_sober + theme(axis.text.x  = element_blank()) +
    coord_cartesian(expand = FALSE, ylim = c(-2, 2), xlim=c(0, 11),clip = "off") + 
  annotate("text", x = 2, y =1.5, label = 'bold(" Evolved\non cherry")',
                     size =3,colour = "#BC3C6D",parse = TRUE) + 
  annotate("text", x = 6, y =1.5, label = 'bold("   Evolved\non cranberry")',
                     size =3,colour = "#FDB424",parse = TRUE) + 
  annotate("text", x = 10, y =1.5, label = 'bold("   Evolved\non strawberry")',
                     size =3,colour = "#3FAA96",parse = TRUE) + 
  ggtitle("Final phenotyping step") 
symp_allop_G29

### Add stroke
DIF_FITNESS_G7 <- symp_allop_g7 + 
  geom_point(aes(alpha = SA, color = interaction(SA,Fruit_s)), position = pd, size =3, stroke =1.5, fill = "white") + 
  scale_alpha_manual(values = c(0, 1)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96", "#7C2748", "#CA8702", "#328677", "#BC3C6D", "#FDB424", "#3FAA96"))  
DIF_FITNESS_G7

DIF_FITNESS_G29 <- symp_allop_G29 + geom_point(aes(alpha = SA, color = interaction(SA,Fruit_s)), position = pd, size =3, stroke =1.5) + 
      scale_alpha_manual(values = c(0, 1)) + 
      scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96", "#7C2748", "#CA8702", "#328677", "#BC3C6D", "#FDB424", "#3FAA96"))  
DIF_FITNESS_G29

plot_legend <- ggplot(data = data_logchange, aes(x = Line, y=logchange, group = Treatment, 
                  color =Fruit_s, shape = Treatment,fill =Fruit_s)) + 
    geom_point(size =3, position = pd) + 
    labs(shape = "Test fruit", color = "Selection fruit") + 
    scale_shape_manual(values = c(21, 22, 24)) + 
    scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96")) + 
    theme_LO_sober + 
    guides(fill = FALSE) + 
    theme(legend.key.size = unit(0.5, "cm"),
        axis.text.x  = element_blank()) +
    coord_cartesian(expand = FALSE, ylim = c(-2, 2), xlim=c(0, 11),clip = "off") 


###########################################################
#############################################ALL PLOT
legend_trade <- lemon::g_legend(plot_legend)
    

SYMP_ALLOP_TOTAL <- cowplot::ggdraw() + 
  cowplot::draw_plot(DIF_FITNESS_G7 + theme(legend.position = 'none', 
                             plot.margin =unit(c(0.5, 3.5, 0.5, 0.5), "cm")),
            x = 0, y = 0.5, width = 1, height = 0.5) + 
  cowplot::draw_plot(DIF_FITNESS_G29 + theme(legend.position = 'none', 
                             plot.margin =unit(c(0.5, 3.5, 0.5, 0.5), "cm")), 
            x = 0, y = 0, width = 1, height = 0.5) + 
  cowplot::draw_plot_label(c("A", "B"),  
                  x = c(0, 0), y = c(1, 0.5), 
                  size = 14) + 
  cowplot::draw_plot(legend_trade, x = 0.93, y = 0.5, width = 0.001, height = 0.001) 

SYMP_ALLOP_TOTAL

cowplot::save_plot(file =here::here("figures", "FIG_SX_HeterogeneityEgg-to-adultviability.pdf"),  
                   SYMP_ALLOP_TOTAL, base_height = 20/cm(1), base_width = 20/cm(1), dpi = 1200)



SYMP_ALLOP_TOTAL <- cowplot::ggdraw() + 
  cowplot::draw_plot(DIF_FITNESS_G7Nb_eggs + theme(legend.position = 'none', 
                             plot.margin =unit(c(0.5, 3.5, 0.5, 0.5), "cm")),
            x = 0, y = 0.5, width = 1, height = 0.5) + 
  cowplot::draw_plot(DIF_FITNESS_G29Nb_eggs + theme(legend.position = 'none', 
                             plot.margin =unit(c(0.5, 3.5, 0.5, 0.5), "cm")), 
            x = 0, y = 0, width = 1, height = 0.5) + 
  cowplot::draw_plot(DIF_FITNESS_G7 + theme(legend.position = 'none', 
                             plot.margin =unit(c(0.5, 3.5, 0.5, 0.5), "cm")),
            x = 0.5, y = 0.5, width = 1, height = 0.5) + 
  cowplot::draw_plot(DIF_FITNESS_G29 + theme(legend.position = 'none', 
                             plot.margin =unit(c(0.5, 3.5, 0.5, 0.5), "cm")), 
            x = 0.5, y = 0, width = 1, height = 0.5) + 
  cowplot::draw_plot_label(c("A", "B", "C", "D"),  
                  x = c(0, 0, 0.5, 0.5), y = c(1, 0.5, 1, 0.5), 
                  size = 14) + 
  cowplot::draw_plot(legend_trade, x = 0.93, y = 0.5, width = 0.001, height = 0.001) 

cowplot::save_plot(file =here::here("figures", "FIG_SX_HeterogeneityNb_Egges_Egg-to-adultviability.pdf"),  
                   SYMP_ALLOP_TOTAL, base_height = 20/cm(1), base_width = 30/cm(1), dpi = 1200)

5 CORRELATED RESPONSES

5.1 Analysis MA and correlation

# ##### Create empty dataframe with slopes and confidence interval for all pairwise
Estimates_pairwise <- data.frame(Variables = rep(rep(c("correlation","cor_CI_inf","cor_CI_sup",
                                                       "slope", "intercept"),2),3),
                                 Generation = rep(c(rep(7, 5), rep(29, 5)),3),
                                 Pairwise = rep(c("Cherry_Cranberry",
                                                     "Cranberry_Strawberry", 
                                                     "Strawberry_Cherry"), each=10),
                                Estimates = NA)

#########################
######################### G7
#########################

####### 
#######  Cherry Cranberry
####### 

## Compute weighted correlation
weightedcor <- sjstats:::weighted_correlation(TEMP_dataG7_CheCran, x=logchange_symp, y=logchange_allop, weights=N_sumsympallop, ci.lvl = 0.95)

## Compute Major Axis regression
model <- lmodel2::lmodel2(logchange_allop ~ logchange_symp,
                            range.y = "interval",range.x = "interval",
                            data = TEMP_dataG7_CheCran[rep(1:nrow(TEMP_dataG7_CheCran), TEMP_dataG7_CheCran$N_sumsympallop),], nperm=0)

## The estimation of the confidence interval with cor.test is wrong
cor.test(x=TEMP_dataG7_CheCran[rep(1:nrow(TEMP_dataG7_CheCran), TEMP_dataG7_CheCran$N_sumsympallop), "logchange_symp"], y=TEMP_dataG7_CheCran[rep(1:nrow(TEMP_dataG7_CheCran), TEMP_dataG7_CheCran$N_sumsympallop), "logchange_allop"])
## 
##  Pearson's product-moment correlation
## 
## data:  TEMP_dataG7_CheCran[rep(1:nrow(TEMP_dataG7_CheCran), TEMP_dataG7_CheCran$N_sumsympallop),  and TEMP_dataG7_CheCran[rep(1:nrow(TEMP_dataG7_CheCran), TEMP_dataG7_CheCran$N_sumsympallop),     "logchange_symp"] and     "logchange_allop"]
## t = 6.3642, df = 136, p-value = 2.79e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3390707 0.5982491
## sample estimates:
##       cor 
## 0.4790334
# Fill in data
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="correlation"] <- weightedcor$estimate
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="cor_CI_inf"] <- weightedcor$ci[1]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="cor_CI_sup"] <- weightedcor$ci[2]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="intercept"] <- model$regression.results[2, 2]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="slope"] <- model$regression.results[2, 3]

####### 
#######  Cranberry Strawberry
####### 

## Compute weighted correlation
weightedcor <- sjstats:::weighted_correlation(TEMP_dataG7_CranStraw, x=logchange_symp, y=logchange_allop, weights=N_sumsympallop, ci.lvl = 0.95)

## Compute Major Axis regression
model <- lmodel2::lmodel2(logchange_allop ~ logchange_symp,
                            range.y = "interval",range.x = "interval",
                            data = TEMP_dataG7_CranStraw[rep(1:nrow(TEMP_dataG7_CranStraw), TEMP_dataG7_CranStraw$N_sumsympallop),], nperm=0)

# Fill in data
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="correlation"] <- weightedcor$estimate
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="cor_CI_inf"] <- weightedcor$ci[1]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="cor_CI_sup"] <- weightedcor$ci[2]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="intercept"] <- model$regression.results[2, 2]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="slope"] <- model$regression.results[2, 3]

####### 
#######  Strawberry  Cherry
####### 

## Compute weighted correlation
weightedcor <- sjstats:::weighted_correlation(TEMP_dataG7_StrawChe, x=logchange_symp, y=logchange_allop, weights=N_sumsympallop, ci.lvl = 0.95)

## Compute Major Axis regression
model <- lmodel2::lmodel2(logchange_allop ~ logchange_symp,
                            range.y = "interval",range.x = "interval",
                            data = TEMP_dataG7_StrawChe[rep(1:nrow(TEMP_dataG7_StrawChe), TEMP_dataG7_StrawChe$N_sumsympallop),], nperm=0)

# Fill in data
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="correlation"] <- weightedcor$estimate
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="cor_CI_inf"] <- weightedcor$ci[1]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="cor_CI_sup"] <- weightedcor$ci[2]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="intercept"] <- model$regression.results[2, 2]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="slope"] <- model$regression.results[2, 3]

#########################
######################### G29
######################### 


####### 
#######  Cherry Cranberry
####### 


## Compute unweighted correlation
unweightedcor <- cor.test(x=TEMP_dataG29_CheCran$logchange_symp, y=TEMP_dataG29_CheCran$logchange_allop)

## Compute weighted correlation
weightedcor <- sjstats:::weighted_correlation(TEMP_dataG29_CheCran, x=logchange_symp, y=logchange_allop, weights=N_sumsympallop, ci.lvl = 0.95)

## Compute Major Axis regression
model <- lmodel2::lmodel2(logchange_allop ~ logchange_symp,
                            range.y = "interval",range.x = "interval",
                            data = TEMP_dataG29_CheCran[rep(1:nrow(TEMP_dataG29_CheCran), TEMP_dataG29_CheCran$N_sumsympallop),], nperm=0)

# Fill in data
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="correlation"] <- unweightedcor$estimate
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="cor_CI_inf"] <- unweightedcor$conf.int[1]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="cor_CI_sup"] <- unweightedcor$conf.int[2]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="intercept"] <- model$regression.results[2, 2]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="slope"] <- model$regression.results[2, 3]

####### 
#######  Cranberry Strawberry
####### 

## Compute unweighted correlation
unweightedcor <- cor.test(x=TEMP_dataG29_CranStraw$logchange_symp, y=TEMP_dataG29_CranStraw$logchange_allop)

## Compute weighted correlation
weightedcor <- sjstats:::weighted_correlation(TEMP_dataG29_CranStraw, x=logchange_symp, y=logchange_allop, weights=N_sumsympallop, ci.lvl = 0.95)

## Compute Major Axis regression
model <- lmodel2::lmodel2(logchange_allop ~ logchange_symp,
                            range.y = "interval",range.x = "interval",
                            data = TEMP_dataG29_CranStraw[rep(1:nrow(TEMP_dataG29_CranStraw), TEMP_dataG29_CranStraw$N_sumsympallop),], nperm=0)

# Fill in data
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="correlation"] <- unweightedcor$estimate
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="cor_CI_inf"] <- unweightedcor$conf.int[1]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="cor_CI_sup"] <- unweightedcor$conf.int[2]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="intercept"] <- model$regression.results[2, 2]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="slope"] <- model$regression.results[2, 3]

####### 
#######  Strawberry  Cherry
####### 

## Compute unweighted correlation
unweightedcor <- cor.test(x=TEMP_dataG29_StrawChe$logchange_symp, y=TEMP_dataG29_StrawChe$logchange_allop)

## Compute weighted correlation
weightedcor <- sjstats:::weighted_correlation(TEMP_dataG29_StrawChe, x=logchange_symp, y=logchange_allop, weights=N_sumsympallop, ci.lvl = 0.95)

## Compute Major Axis regression
model <- lmodel2::lmodel2(logchange_allop ~ logchange_symp,
                            range.y = "interval",range.x = "interval",
                            data = TEMP_dataG29_StrawChe[rep(1:nrow(TEMP_dataG29_StrawChe), TEMP_dataG29_StrawChe$N_sumsympallop),], nperm=0)

# Fill in data
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="correlation"] <- unweightedcor$estimate
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="cor_CI_inf"] <- unweightedcor$conf.int[1]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="cor_CI_sup"] <- unweightedcor$conf.int[2]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="intercept"] <- model$regression.results[2, 2]
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="slope"] <- model$regression.results[2, 3]



Estimates_pairwise
##      Variables Generation             Pairwise    Estimates
## 1  correlation          7     Cherry_Cranberry  0.479033353
## 2   cor_CI_inf          7     Cherry_Cranberry  0.001431407
## 3   cor_CI_sup          7     Cherry_Cranberry  0.956635299
## 4        slope          7     Cherry_Cranberry  4.160858972
## 5    intercept          7     Cherry_Cranberry -0.356412153
## 6  correlation         29     Cherry_Cranberry -0.567557085
## 7   cor_CI_inf         29     Cherry_Cranberry -0.908773300
## 8   cor_CI_sup         29     Cherry_Cranberry  0.228504360
## 9        slope         29     Cherry_Cranberry -1.427127889
## 10   intercept         29     Cherry_Cranberry  0.268405728
## 11 correlation          7 Cranberry_Strawberry  0.535360085
## 12  cor_CI_inf          7 Cranberry_Strawberry  0.105487696
## 13  cor_CI_sup          7 Cranberry_Strawberry  0.965232474
## 14       slope          7 Cranberry_Strawberry  1.070761415
## 15   intercept          7 Cranberry_Strawberry -0.049525449
## 16 correlation         29 Cranberry_Strawberry  0.593494172
## 17  cor_CI_inf         29 Cranberry_Strawberry -0.191100790
## 18  cor_CI_sup         29 Cranberry_Strawberry  0.915350061
## 19       slope         29 Cranberry_Strawberry  1.613908561
## 20   intercept         29 Cranberry_Strawberry -0.582936853
## 21 correlation          7    Strawberry_Cherry  0.616674051
## 22  cor_CI_inf          7    Strawberry_Cherry  0.188354905
## 23  cor_CI_sup          7    Strawberry_Cherry  1.044993196
## 24       slope          7    Strawberry_Cherry  0.829196290
## 25   intercept          7    Strawberry_Cherry  0.177499018
## 26 correlation         29    Strawberry_Cherry  0.676000585
## 27  cor_CI_inf         29    Strawberry_Cherry -0.300322260
## 28  cor_CI_sup         29    Strawberry_Cherry  0.960575095
## 29       slope         29    Strawberry_Cherry  2.537639762
## 30   intercept         29    Strawberry_Cherry -0.764708220

5.2 Compute confidence interval of the difference

computeCIcordifG29G7(pair="Cherry_Cranberry")
##        dif      lower     higher 
## -1.0465904 -1.9126980 -0.1804829
computeCIcordifG29G7(pair="Cranberry_Strawberry")
##         dif       lower      higher 
##  0.05813409 -0.78991100  0.90617917
computeCIcordifG29G7(pair="Strawberry_Cherry")
##         dif       lower      higher 
##  0.05932653 -0.95762427  1.07627734
# ##### Create empty dataframe with slopes and confidence interval for all pairwise
Estimates_pairwise <- data.frame(Variables = rep(rep(c("correlation","cor_CI_inf","cor_CI_sup",
                                                       "slope", "intercept"),2),3),
                                 Generation = rep(c(rep(7, 5), rep(29, 5)),3),
                                 Pairwise = rep(c("Cherry_Cranberry",
                                                     "Cranberry_Strawberry", 
                                                     "Strawberry_Cherry"), each=10),
                                Estimates = NA ,
                                EstimatesWeighed = NA ,
                                EstimatesUnweighed = NA )

# Calculate sd in case of NA value 
mean_sd_G7 <- sqrt(mean(data_logchange$sd_logchange[data_logchange$Generation=="7"]^2, na.rm = TRUE))
mean_sd_G29 <- sqrt(mean(data_logchange$sd_logchange[data_logchange$Generation=="29"]^2, na.rm = TRUE))

## Check whether the squared standard errors are of the same magnitude across fruit media (answer:yes)
tapply(data_logchange$sd_logchange^2, list(data_logchange$Treatment, data_logchange$Generation), mean, na.rm = TRUE)

## Impute directly the missing standard errors
data_logchange$se_logchange <- data_logchange$sd_logchange
data_logchange$se_logchange[data_logchange$Generation=="7"&is.na(data_logchange$se_logchange)] <- rep(mean_sd_G7, sum(data_logchange$Generation=="7"&is.na(data_logchange$se_logchange)))


## Nbr sample 
nb_sim = 5000


#########################
######################### G7
#########################

####### 
#######  Cherry Cranberry
####### 


# MA axis regression with bootstrap to estimate slope, intercept and 
sim <- data.frame(t(sapply(1:nb_sim, bootstrap_lmodel2, data_pairwise = TEMP_dataG7_CheCran, sd_NA = mean_sd_G7)))

# Fill out data with quantile of simulations
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="correlation"] <- mean(sim$correlation.cor, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="cor_CI_inf"] <- quantile(sim$correlation.cor, probs = 0.025, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="cor_CI_sup"] <- quantile(sim$correlation.cor, probs = 0.975, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="intercept"] <- mean(sim$intercept, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="slope"] <- mean(sim$slope, na.rm = TRUE)

## Add other correlations
Estimates_pairwise$EstimatesWeighed[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="correlation"] <- weighted_correlation(TEMP_dataG7_CheCran, x=logchange_symp, y=logchange_allop, weights=N_sumsympallop, ci.lvl = 0.95)$estimate

Estimates_pairwise$EstimatesWeighed[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="cor_CI_inf"] <- weighted_correlation(TEMP_dataG7_CheCran, x=logchange_symp, y=logchange_allop, weights=N_sumsympallop, ci.lvl = 0.95)$ci[1]

Estimates_pairwise$EstimatesWeighed[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="cor_CI_sup"] <- weighted_correlation(TEMP_dataG7_CheCran, x=logchange_symp, y=logchange_allop, weights=N_sumsympallop, ci.lvl = 0.95)$ci[2]

Estimates_pairwise$EstimatesUnweighed[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="correlation"] <- cor.test(x=TEMP_dataG7_CheCran$logchange_symp, y=TEMP_dataG7_CheCran$logchange_allop)$estimate

Estimates_pairwise$EstimatesUnweighed[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="cor_CI_inf"] <- cor.test(x=TEMP_dataG7_CheCran$logchange_symp, y=TEMP_dataG7_CheCran$logchange_allop)$conf.int[1]

Estimates_pairwise$EstimatesUnweighed[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="cor_CI_sup"] <- cor.test(x=TEMP_dataG7_CheCran$logchange_symp, y=TEMP_dataG7_CheCran$logchange_allop)$conf.int[2]



#Save sim G7 to calculate proportion later: 
sim_CheCranG7<-sim

####### 
#######  Cranberry Strawberry
####### 

# MA axis regression with bootstrap to estimate slope, intercept and 
sim <- data.frame(t(sapply(1:nb_sim, bootstrap_lmodel2, data_pairwise = TEMP_dataG7_CranStraw, sd_NA = mean_sd_G7)))

# Fill out data with quantile of simulations
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="correlation"] <- quantile(sim$cor, probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="cor_CI_inf"] <- quantile(sim$cor, probs = 0.025, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="cor_CI_sup"] <- quantile(sim$cor, probs = 0.975, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="intercept"] <- quantile(sim$intercept,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="slope"] <- quantile(sim$slope,probs = 0.5, na.rm = TRUE)


## Save sim G7 to calculate correlation difference later: 
sim_CranStrawG7<-sim



####### 
#######  Strawberry  Cherry
####### 



# MA axis regression with bootstrap to estimate slope, intercept and 
sim <- data.frame(t(sapply(1:nb_sim, bootstrap_lmodel2, data_pairwise = TEMP_dataG7_StrawChe, sd_NA = mean_sd_G7)))

# Fill out data with quantile of simulations
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="correlation"] <- quantile(sim$cor, probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="cor_CI_inf"] <- quantile(sim$cor, probs = 0.025, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="cor_CI_sup"] <- quantile(sim$cor, probs = 0.975, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="intercept"] <- quantile(sim$intercept,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==7&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="slope"] <- quantile(sim$slope,probs = 0.5, na.rm = TRUE)


## Save sim G7 to calculate correlation difference later: 
sim_StrawCheG7<-sim



#########################
######################### G29
######################### 


####### 
#######  Cherry Cranberry
####### 

# MA axis regression with bootstrap to estimate slope, intercept and 
sim <- data.frame(t(sapply(1:nb_sim, bootstrap_lmodel2, data_pairwise = TEMP_dataG29_CheCran, sd_NA = mean_sd_G29)))

# Fill out data with quantile of simulations
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="correlation"] <- quantile(sim$cor, probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="cor_CI_inf"] <- quantile(sim$cor, probs = 0.025, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="cor_CI_sup"] <- quantile(sim$cor, probs = 0.975, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="intercept"] <- quantile(sim$intercept, probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="slope"] <- quantile(sim$slope, probs = 0.5, na.rm = TRUE)


## Save sim G29 to calculate correlation difference later: 
sim_CheCranG29<-sim

####### 
#######  Cranberry Strawberry
####### 

# MA axis regression with bootstrap to estimate slope, intercept and 
sim <- data.frame(t(sapply(1:nb_sim, bootstrap_lmodel2, data_pairwise = TEMP_dataG29_CranStraw, sd_NA = mean_sd_G29)))

# Fill out data with quantile of simulations
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="correlation"] <- quantile(sim$cor, probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="cor_CI_inf"] <- quantile(sim$cor, probs = 0.025, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="cor_CI_sup"] <- quantile(sim$cor, probs = 0.975, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="intercept"] <- quantile(sim$intercept,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cranberry_Strawberry"&
                              Estimates_pairwise$Variables=="slope"] <- quantile(sim$slope,probs = 0.5, na.rm = TRUE)

## Save sim G29 to calculate correlation difference later: 
sim_CranStrawG29<-sim

####### 
#######  Strawberry  Cherry
####### 
rm(sim)
# MA axis regression with bootstrap to estimate slope, intercept and 
sim <- data.frame(t(sapply(1:nb_sim, bootstrap_lmodel2, data_pairwise = TEMP_dataG29_StrawChe, sd_NA = mean_sd_G29)))

# Fill out data with quantile of simulations
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="correlation"] <- quantile(sim$cor, probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="cor_CI_inf"] <- quantile(sim$cor, probs = 0.025, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="cor_CI_sup"] <- quantile(sim$cor, probs = 0.975, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="intercept"] <- quantile(sim$intercept,probs = 0.5, na.rm = TRUE)
Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Strawberry_Cherry"&
                              Estimates_pairwise$Variables=="slope"] <- quantile(sim$slope,probs = 0.5, na.rm = TRUE)



Estimates_pairwise

## Save sim G29 to calculate correlation difference later: 
sim_StrawCheG29<-sim

## Compute the confidence interval of the difference between correlations in G7 and G29
difcorCheCran <- sim_CheCranG29$correlation.cor-sim_CheCranG7$correlation.cor
quantile(difcorCheCran, probs = c(0.025, 0.975), na.rm=TRUE)

difcorCranStraw <- sim_CranStrawG29$correlation.cor-sim_CranStrawG7$correlation.cor
quantile(difcorCranStraw, probs = c(0.025, 0.975), na.rm=TRUE)

difcorStrawChe <- sim_StrawCheG29$correlation.cor-sim_StrawCheG7$correlation.cor
quantile(difcorStrawChe, probs = c(0.025, 0.975), na.rm=TRUE)

#Cherry Cranbery 
#Compute proportion of simulations correlations G7 are lower than correlation estimated on G29 
cor_G29_CheCran<-Estimates_pairwise$Estimates[Estimates_pairwise$Generation==29&
                              Estimates_pairwise$Pairwise=="Cherry_Cranberry"&
                              Estimates_pairwise$Variables=="correlation"]
#Estimates of correlation during final phenotyping for Cherry Cranberry 
p_val_decreaseCheCran<-length(sim_CheCranG7$correlation.cor[sim_CheCranG7$correlation.cor<cor_G29_CheCran])/nb_sim


p_val_decreaseCheCran

5.3 Plot

ymin = -50
ymax = 50


## Function equation R
eq_r <- function(gen = 7, pair = "Cherry_Cranberry") {
  if(gen == 7 &pair == "Cherry_Cranberry"){
            as.character(
      as.expression(
        substitute(~~italic(r)^2~"="~r2~"["~inf~","~sup~"]",
                  list(r2 = format(Estimates_pairwise$Estimates[Estimates_pairwise$Generation == gen &
  Estimates_pairwise$Pairwise == pair &
  Estimates_pairwise$Variables == "correlation"], digits = 2), 
                  inf = format(Estimates_pairwise$Estimates[Estimates_pairwise$Generation == gen &
  Estimates_pairwise$Pairwise == pair &
  Estimates_pairwise$Variables == "cor_CI_inf"], digits = 1, nsmall=2), 
                  sup = format(Estimates_pairwise$Estimates[Estimates_pairwise$Generation == gen &
  Estimates_pairwise$Pairwise == pair &
  Estimates_pairwise$Variables == "cor_CI_sup"], digits = 2)))
      )
    )
      

  }else{
      as.character(
    as.expression(
      substitute(~~italic(r)^2~"="~r2~"["~inf~","~sup~"]",
                list(r2 = format(Estimates_pairwise$Estimates[Estimates_pairwise$Generation == gen &
Estimates_pairwise$Pairwise == pair &
Estimates_pairwise$Variables == "correlation"], digits = 2), 
                inf = format(Estimates_pairwise$Estimates[Estimates_pairwise$Generation == gen &
Estimates_pairwise$Pairwise == pair &
Estimates_pairwise$Variables == "cor_CI_inf"], digits = 2, nsmall=2), 
                sup = format(Estimates_pairwise$Estimates[Estimates_pairwise$Generation == gen &
Estimates_pairwise$Pairwise == pair &
Estimates_pairwise$Variables == "cor_CI_sup"], digits = 2)))
    )
  )
  }

}



##################################################
##################   G7
# Limits
ymin_CheCranG7=min(min(TEMP_dataG7_CheCran$logchange_allop-1.96*TEMP_dataG7_CheCran$sd_allop, na.rm= TRUE),
                min(TEMP_dataG7_CheCran$logchange_symp-1.96*TEMP_dataG7_CheCran$sd_symp, na.rm= TRUE))
ymax_CheCranG7=max(max(TEMP_dataG7_CheCran$logchange_allop + 1.96*TEMP_dataG7_CheCran$sd_allop, na.rm= TRUE),
                max(TEMP_dataG7_CheCran$logchange_symp + 1.96*TEMP_dataG7_CheCran$sd_symp, na.rm= TRUE))
lim_text<-ymin_CheCranG7+0.99*(ymax_CheCranG7-ymin_CheCranG7)



# Plot
CheCran_G7 <- ggplot(data = TEMP_dataG7_CheCran) + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "slope"],
              colour = "black", size = 0.75) + 
  geom_text(x = -0.825, y = lim_text, label = eq_r(gen = 7, pair = "Cherry_Cranberry"), parse = TRUE, color="black", size = 3.5) +
  geom_errorbar(aes(x = logchange_symp, ymin = logchange_allop-(1.96*sd_allop),
                    ymax = logchange_allop + (1.96*sd_allop),
                    color = Symp, linetype = Line_type),
                width= 0.02, size = 0.3, alpha = 0.8) + 
  geom_errorbarh(aes(y = logchange_allop, xmin = logchange_symp-1.96*sd_symp, xmax = logchange_symp + 1.96*sd_symp, 
                     color = Symp, linetype = Line_type),
                 height = 0.02, size = 0.3, alpha = 0.8) + 
  geom_point(aes(x = logchange_symp, y = logchange_allop,  color = Symp,fill = Symp, shape = Allop),
             size = 3, fill = "white", stroke =1.2) + 
  xlab("Fitness change in\nselective environment")  + 
  ylab("Fitness change in\nalternative environment")  + 
  ggtitle("Cherry vs. Cranberry") + 
  theme_LO_sober + 
  labs(shape = "Test fruit", color = "Evolution fruit") + 
  scale_shape_manual(values = c(21, 22)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424"))  + 
  coord_cartesian(ylim = c(ymin_CheCranG7, ymax_CheCranG7), 
                  xlim = c(ymin_CheCranG7, ymax_CheCranG7)) 
  
 CheCran_G7

# Limits
ymin_CranStrawG7=min(min(TEMP_dataG7_CranStraw$logchange_allop-1.96*TEMP_dataG7_CranStraw$sd_allop, na.rm= TRUE),
                min(TEMP_dataG7_CranStraw$logchange_symp-1.96*TEMP_dataG7_CranStraw$sd_symp, na.rm= TRUE))
ymax_CranStrawG7=max(max(TEMP_dataG7_CranStraw$logchange_allop + 1.96*TEMP_dataG7_CranStraw$sd_allop, na.rm= TRUE),
                max(TEMP_dataG7_CranStraw$logchange_symp + 1.96*TEMP_dataG7_CranStraw$sd_symp, na.rm= TRUE))
lim_text<-ymin_CranStrawG7+0.99*(ymax_CranStrawG7-ymin_CranStrawG7)


CranStraw_G7 <-  ggplot(data = TEMP_dataG7_CranStraw) + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry" & 
                               Estimates_pairwise$Variables == "slope"],
              colour = "black", size = 0.75) + 
  geom_errorbar(aes(x =logchange_symp, ymin = logchange_allop-(1.96*sd_allop), 
                    ymax = logchange_allop + (1.96*sd_allop),
                    color = Symp, linetype = Line_type),
                  width= 0.02, size = 0.3, alpha = 0.8) + 
  geom_errorbarh(aes(y = logchange_allop, xmin = logchange_symp-1.96*sd_symp, xmax = logchange_symp + 1.96*sd_symp, 
                 color = Symp, linetype = Line_type),
                 height = 0.02, size = 0.3, alpha = 0.8) + 
  geom_point(aes(x =logchange_symp, y = logchange_allop,  color = Symp,fill = Symp, shape = Allop),
                   size =3, fill = "white", stroke =1.2) + 
  geom_text(x = -0.45, y = lim_text, label = eq_r(gen = 7, pair = "Cranberry_Strawberry"), parse = TRUE, color="black", size = 3.5) +
  xlab("Fitness change in\nselective environment")  + 
  ylab("Fitness change in\nalternative environment")  + 
     ggtitle("Cranberry vs. Strawberry") + 
  coord_cartesian(ylim = c(ymin_CranStrawG7, ymax_CranStrawG7), 
                  xlim = c(ymin_CranStrawG7, ymax_CranStrawG7)) + 
   labs(shape = "Test fruit", color = "Selection fruit") + 
   scale_shape_manual(values = c(22, 24)) + 
   scale_color_manual(values = c("#FDB424", "#3FAA96"))  + 
  theme_LO_sober
 CranStraw_G7

# Limits
ymin_StrawCheG7=min(min(TEMP_dataG7_StrawChe$logchange_allop-1.96*TEMP_dataG7_StrawChe$sd_allop, na.rm= TRUE),
                min(TEMP_dataG7_StrawChe$logchange_symp-1.96*TEMP_dataG7_StrawChe$sd_symp, na.rm= TRUE))
ymax_StrawCheG7=max(max(TEMP_dataG7_StrawChe$logchange_allop + 1.96*TEMP_dataG7_StrawChe$sd_allop, na.rm= TRUE),
                max(TEMP_dataG7_StrawChe$logchange_symp + 1.96*TEMP_dataG7_StrawChe$sd_symp, na.rm= TRUE))
lim_text<-ymin_StrawCheG7+0.99*(ymax_StrawCheG7-ymin_StrawCheG7)


StrawChe_G7 <-  ggplot(data = TEMP_dataG7_StrawChe) + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation == 7 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "slope"],
              colour = "black", size = 0.75) + 
  geom_errorbar(aes(x =logchange_symp, ymin = logchange_allop-(1.96*sd_allop), ymax = logchange_allop + (1.96*sd_allop),
                    color = Symp, linetype = Line_type),
                  width= 0.02, size = 0.3, alpha = 0.8) + 
  geom_errorbarh(aes(y = logchange_allop, xmin = logchange_symp-1.96*sd_symp, xmax = logchange_symp + 1.96*sd_symp, 
                 color = Symp, linetype = Line_type),
                 height = 0.02, size = 0.3, alpha = 0.8) + 
  geom_point(aes(x = logchange_symp, y = logchange_allop,  color = Symp,fill = Symp, shape = Allop),
                 size =3, fill = "white", stroke =1.2) + 
  geom_text(x = -0.83, y = lim_text, label = eq_r(gen = 7, pair = "Strawberry_Cherry"), parse = TRUE, color="black", size = 3.5) +
  xlab("Fitness change in\nselective environment")  + 
  ylab("Fitness change in\nalternative environment")  + 
  ggtitle("Strawberry vs. Cherry") + 
  labs(shape = "Test fruit", color = "Selection fruit") + 
  scale_shape_manual(values = c(21, 24)) + 
  scale_color_manual(values = c("#BC3C6D", "#3FAA96"))  + 
  coord_cartesian(ylim = c(ymin_StrawCheG7, ymax_StrawCheG7), 
                  xlim = c(ymin_StrawCheG7, ymax_StrawCheG7)) + 
  theme_LO_sober
 StrawChe_G7

#####################################################################################
##################################      G29        ##################################
#####################################################################################
TEMP_dataG29_CheCran
##    Line Treatment   Fruit_s Generation N_allop   logchange sd_logchange SA
## 15  CEA Cranberry    Cherry         29      30 -0.21774071   0.13281291  0
## 17  CEB Cranberry    Cherry         29      30 -0.59200204   0.09634940  0
## 21  CEC Cranberry    Cherry         29      30 -0.20670829   0.12591289  0
## 37  CRA    Cherry Cranberry         29      30  0.32537794   0.10025669  0
## 41  CRB    Cherry Cranberry         29      30  0.16929862   0.12699030  0
## 45  CRC    Cherry Cranberry         29      30  0.08474123   0.12334598  0
## 47  CRD    Cherry Cranberry         29      30  0.54736475   0.09718018  0
## 49  CRE    Cherry Cranberry         29      30 -0.10496044   0.11556266  0
##         Symp     Allop Line_type logchange_allop   sd_allop logchange_symp
## 15    Cherry Cranberry    dashed     -0.21774071 0.13281291     0.14767000
## 17    Cherry Cranberry    dashed     -0.59200204 0.09634940     0.25038029
## 21    Cherry Cranberry    dashed     -0.20670829 0.12591289     0.35854062
## 37 Cranberry    Cherry    dashed      0.32537794 0.10025669    -0.14885350
## 41 Cranberry    Cherry    dashed      0.16929862 0.12699030     0.29719689
## 45 Cranberry    Cherry    dashed      0.08474123 0.12334598     0.02703717
## 47 Cranberry    Cherry    dashed      0.54736475 0.09718018    -0.15627744
## 49 Cranberry    Cherry    dashed     -0.10496044 0.11556266     0.72513486
##      sd_symp N_symp N_sumsympallop
## 15 0.1574284     30             60
## 17 0.1009864     30             60
## 21 0.1169300     30             60
## 37 0.1193552     30             60
## 41 0.1531197     30             60
## 45 0.1156381     30             60
## 47 0.1156695     30             60
## 49 0.1275430     30             60
TEMP_dataG29_CranStraw
##    Line  Treatment    Fruit_s Generation N_allop   logchange sd_logchange SA
## 39  CRA Strawberry  Cranberry         29      30 -0.29731422    0.1098141  0
## 42  CRB Strawberry  Cranberry         29      30 -0.24539416    0.1387043  0
## 43  CRC Strawberry  Cranberry         29      30 -0.28950778    0.2012393  0
## 46  CRD Strawberry  Cranberry         29      30 -1.02555272    0.1902082  0
## 51  CRE Strawberry  Cranberry         29      30 -0.16810249    0.1398057  0
## 67  FRA  Cranberry Strawberry         29      30 -0.25158762    0.1712755  0
## 72  FRB  Cranberry Strawberry         29      30  0.04303751    0.1441119  0
## 73  FRC  Cranberry Strawberry         29      30  0.52109773    0.1254820  0
##          Symp      Allop Line_type logchange_allop  sd_allop logchange_symp
## 39  Cranberry Strawberry    dashed     -0.29731422 0.1098141    -0.14885350
## 42  Cranberry Strawberry    dashed     -0.24539416 0.1387043     0.29719689
## 43  Cranberry Strawberry    dashed     -0.28950778 0.2012393     0.02703717
## 46  Cranberry Strawberry    dashed     -1.02555272 0.1902082    -0.15627744
## 51  Cranberry Strawberry    dashed     -0.16810249 0.1398057     0.72513486
## 67 Strawberry  Cranberry    dashed     -0.25158762 0.1712755     0.15709032
## 72 Strawberry  Cranberry    dashed      0.04303751 0.1441119     0.56023100
## 73 Strawberry  Cranberry    dashed      0.52109773 0.1254820     0.36640738
##      sd_symp N_symp N_sumsympallop
## 39 0.1193552     30             60
## 42 0.1531197     30             60
## 43 0.1156381     30             60
## 46 0.1156695     30             60
## 51 0.1275430     30             60
## 67 0.1411560     30             60
## 72 0.1284283     30             60
## 73 0.1201170     30             60
TEMP_dataG29_StrawChe
##    Line  Treatment    Fruit_s Generation N_allop   logchange sd_logchange SA
## 14  CEA Strawberry     Cherry         29      30  0.01971359    0.1475176  0
## 16  CEB Strawberry     Cherry         29      30 -0.15850680    0.0983504  0
## 19  CEC Strawberry     Cherry         29      30  0.21131464    0.1237821  0
## 68  FRA     Cherry Strawberry         29      30 -0.50599095    0.1378652  0
## 71  FRB     Cherry Strawberry         29      30  0.19377147    0.1094819  0
## 75  FRC     Cherry Strawberry         29      30  0.32151693    0.1140574  0
##          Symp      Allop Line_type logchange_allop  sd_allop logchange_symp
## 14     Cherry Strawberry    dashed      0.01971359 0.1475176      0.1476700
## 16     Cherry Strawberry    dashed     -0.15850680 0.0983504      0.2503803
## 19     Cherry Strawberry    dashed      0.21131464 0.1237821      0.3585406
## 68 Strawberry     Cherry    dashed     -0.50599095 0.1378652      0.1570903
## 71 Strawberry     Cherry    dashed      0.19377147 0.1094819      0.5602310
## 75 Strawberry     Cherry    dashed      0.32151693 0.1140574      0.3664074
##      sd_symp N_symp N_sumsympallop
## 14 0.1574284     30             60
## 16 0.1009864     30             60
## 19 0.1169300     30             60
## 68 0.1411560     30             60
## 71 0.1284283     30             60
## 75 0.1201170     30             60
# Limits
ymin_CheCranG29=min(min(TEMP_dataG29_CheCran$logchange_allop-1.96*TEMP_dataG29_CheCran$sd_allop, na.rm= TRUE),
                min(TEMP_dataG29_CheCran$logchange_symp-1.96*TEMP_dataG29_CheCran$sd_symp, na.rm= TRUE))
ymax_CheCranG29=max(max(TEMP_dataG29_CheCran$logchange_allop + 1.96*TEMP_dataG29_CheCran$sd_allop, na.rm= TRUE),
                max(TEMP_dataG29_CheCran$logchange_symp + 1.96*TEMP_dataG29_CheCran$sd_symp, na.rm= TRUE))
lim_text<-ymin_CheCranG29+0.99*(ymax_CheCranG29-ymin_CheCranG29)


CheCran_G29 <- ggplot(data = TEMP_dataG29_CheCran) + 
    geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cherry_Cranberry" & 
                               Estimates_pairwise$Variables == "slope"],
              colour = "black", size = 0.75) + 
  geom_errorbar(aes(x =logchange_symp, ymin = logchange_allop-(1.96*sd_allop), 
                    ymax = logchange_allop + (1.96*sd_allop),
                    color = Symp),
                  width= 0.02, size = 0.3, alpha = 0.8) + 
  geom_errorbarh(aes(y = logchange_allop, xmin = logchange_symp-1.96*sd_symp, xmax = logchange_symp + 1.96*sd_symp, 
                 color = Symp),
                 height = 0.02, size = 0.3, alpha = 0.8) + 
  geom_point(aes(x =logchange_symp, y = logchange_allop,  color = Symp,fill = Symp, shape = Allop),
                 size =3, fill = "white", stroke =1.2) + 
  geom_text(x = 0.33, y = lim_text, label = eq_r(gen = 29, pair = "Cherry_Cranberry"), parse = TRUE, color="black", size = 3.5) +
  xlab("Fitness change in\nselective environment")  + 
  ylab("Fitness change in\nalternative environment")  + 
  ggtitle("Cherry vs. Cranberry") + 
  labs(shape = "Test fruit", color = "Evolution fruit") + 
  scale_shape_manual(values = c(16, 15)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424"))  + 
  coord_cartesian(ylim = c(ymin_CheCranG29, ymax_CheCranG29), 
                  xlim = c(ymin_CheCranG29, ymax_CheCranG29)) + 
  theme_LO_sober
 CheCran_G29

# Limits
ymin_CranStrawG29=min(min(TEMP_dataG29_CranStraw$logchange_allop-1.96*TEMP_dataG29_CranStraw$sd_allop, na.rm= TRUE),
                min(TEMP_dataG29_CranStraw$logchange_symp-1.96*TEMP_dataG29_CranStraw$sd_symp, na.rm= TRUE))
ymax_CranStrawG29=max(max(TEMP_dataG29_CranStraw$logchange_allop + 1.96*TEMP_dataG29_CranStraw$sd_allop, na.rm= TRUE),
                max(TEMP_dataG29_CranStraw$logchange_symp + 1.96*TEMP_dataG29_CranStraw$sd_symp, na.rm= TRUE))
lim_text<-ymin_CranStrawG29+0.99*(ymax_CranStrawG29-ymin_CranStrawG29)

#Plot
 CranStraw_G29 <- ggplot(data = TEMP_dataG29_CranStraw) + 
  geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Cranberry_Strawberry" & 
                               Estimates_pairwise$Variables == "slope"],
              colour = "black", size = 0.75) + 
   geom_errorbar(aes(x =logchange_symp, ymin = logchange_allop-(1.96*sd_allop), 
                    ymax = logchange_allop + (1.96*sd_allop),
                    color = Symp),
                  width= 0.02, size = 0.3, alpha = 0.8) + 
   geom_errorbarh(aes(y = logchange_allop, xmin = logchange_symp-1.96*sd_symp, xmax = logchange_symp + 1.96*sd_symp, 
                 color = Symp),
                 height = 0.02, size = 0.3, alpha = 0.8) + 
   geom_point(aes(x =logchange_symp, y = logchange_allop,  color = Symp,fill = Symp, shape = Allop),
                 size =3, fill = "white", stroke =1.2) + 
   geom_text(x = -0.45, y = lim_text, label = eq_r(gen = 29, pair = "Cranberry_Strawberry"), parse = TRUE, color="black", size = 3.5) +
   xlab("Fitness change in\nselective environment")  + 
   ylab("Fitness change in\nalternative environment")  + 
   ggtitle("Cranberry vs. Strawberry") + 
   coord_cartesian(ylim = c(ymin_CranStrawG29, ymax_CranStrawG29), 
                  xlim = c(ymin_CranStrawG29, ymax_CranStrawG29)) + 
   labs(shape = "Test fruit", color = "Selection fruit") + 
   scale_shape_manual(values = c(15, 17)) + 
   scale_color_manual(values = c("#FDB424", "#3FAA96"))  + 
   theme_LO_sober
 CranStraw_G29

# Limits
ymin_StrawCheG29=min(min(TEMP_dataG29_StrawChe$logchange_allop-1.96*TEMP_dataG29_StrawChe$sd_allop, na.rm= TRUE),
                min(TEMP_dataG29_StrawChe$logchange_symp-1.96*TEMP_dataG29_StrawChe$sd_symp, na.rm= TRUE))
ymax_StrawCheG29=max(max(TEMP_dataG29_StrawChe$logchange_allop + 1.96*TEMP_dataG29_StrawChe$sd_allop, na.rm= TRUE),
                max(TEMP_dataG29_StrawChe$logchange_symp + 1.96*TEMP_dataG29_StrawChe$sd_symp, na.rm= TRUE))
lim_text<-ymin_StrawCheG29+0.99*(ymax_StrawCheG29-ymin_StrawCheG29)


StrawChe_G29 <- ggplot(data = TEMP_dataG29_StrawChe) + 
 geom_abline(intercept = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "intercept"],
               slope = Estimates_pairwise$Estimates[Estimates_pairwise$Generation ==29 & 
                               Estimates_pairwise$Pairwise == "Strawberry_Cherry" & 
                               Estimates_pairwise$Variables == "slope"],
              colour = "black", size = 0.75) + 
  geom_errorbar(aes(x =logchange_symp, ymin = logchange_allop-(1.96*sd_allop), 
                    ymax = logchange_allop + (1.96*sd_allop),
                    color = Symp),
                  width= 0.02, size = 0.2, alpha =1) + 
  geom_errorbarh(aes(y = logchange_allop, xmin = logchange_symp-1.96*sd_symp, xmax = logchange_symp + 1.96*sd_symp, 
                 color = Symp),
                 height = 0.02, size = 0.2, alpha =1) + 
  geom_point(aes(x =logchange_symp, y = logchange_allop,  color = Symp,fill = Symp, shape = Allop),
                 size =3, fill = "white", stroke =1.2) + 
  geom_text(x = -0.2, y = lim_text, label = eq_r(gen = 29, pair = "Strawberry_Cherry"), parse = TRUE, color="black", size = 3.5) +
  coord_cartesian(ylim = c(ymin_StrawCheG29, ymax_StrawCheG29), 
                  xlim = c(ymin_StrawCheG29, ymax_StrawCheG29)) + 
    xlab("Fitness change in\nselective environment")  + 
    ylab("Fitness change in\nalternative environment")  + 
  ggtitle("Strawberry vs. Cherry") + 
  labs(shape = "Test fruit", color = "Selection fruit") + 
  scale_shape_manual(values = c(16, 17)) + 
  scale_color_manual(values = c("#BC3C6D", "#3FAA96"))  + 
  theme_LO_sober
 StrawChe_G29

legend_tradevide <-  ggplot(data = data_logchange[data_logchange$Generation == "29",],
                          aes(x =logchange, y = logchange,  color = Symp,fill = Symp, shape = Allop)) + 
  geom_point(size =2.5, fill = "white") + 
  labs(shape = "Test fruit", color = "Selection fruit") + 
  scale_shape_manual(values = c(16, 15, 17)) + 
  scale_color_manual(values = c("#BC3C6D", "#FDB424", "#3FAA96"))  + 
  theme_LO_sober

legend_trade <- lemon::g_legend(legend_tradevide)

 
 

Slopeestimates_logchange_pairwise_errorbar <- cowplot::ggdraw() + 
  cowplot::draw_plot(CheCran_G7 + theme(legend.position = "none"), 
            x = 0.01, y = 0.5, width = 0.22, height = 0.45) + 
  cowplot::draw_plot(CranStraw_G7 + theme(legend.position = "none"), 
            x = 0.31, y = 0.5, width = 0.22, height = 0.45) + 
  cowplot::draw_plot(StrawChe_G7 + theme(legend.position = "none"), 
            x = 0.61, y = 0.5, width = 0.22, height = 0.45) + 
  cowplot::draw_plot(legend_trade, x = 0.85, y = 0.5, width = 0.1, height = 0.1) + 
  cowplot::draw_plot(CheCran_G29 + theme(legend.position = "none"), 
            x = 0.01, y = 0, width = 0.22, height = 0.45) + 
  cowplot::draw_plot(CranStraw_G29 + theme(legend.position = "none"), 
            x = 0.31, y = 0, width = 0.22, height = 0.45) + 
  cowplot::draw_plot(StrawChe_G29 + theme(legend.position = "none"), 
            x = 0.61, y = 0, width = 0.22, height = 0.45) + 
  cowplot::draw_plot_label(c("Intermediate phenotyping step", "A", "B", "C", " ",
                    "Final phenotyping step", "D", "E", "F", " "),  
                  x = c(0.30, 0.01, 0.30, 0.61, 0.92, 0.250, 0.01, 0.30, 0.61, 0.92), 
                  y = c(0.99, 0.95, 0.95, 0.95, 0.95, 0.49, 0.45, 0.45, 0.45, 0.45), 
                  hjust = c(-0.25, -0.25, -0.25, -0.25, -0.25, -0.75, -0.75, -0.75, -0.75, -0.75), 
                  vjust = c(1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5),
                  size = 16) 
 Slopeestimates_logchange_pairwise_errorbar

#  
# #  
# 
cowplot::save_plot(file =here::here("figures", "FIG_Correlation.pdf"),
                   Slopeestimates_logchange_pairwise_errorbar,
                   base_height = 18/cm(1), base_width = 34/cm(1), dpi = 610)

6 SA AND LOCAL ADAPTATION

6.1 Analysis fitness change

# Add vector of weight 
data_logchange$Inv_sum<-NA
data_logchange$Vector_sample<-NA
data_logchange[data_logchange$Generation=="7",]$Inv_sum <- 1 / (data_logchange[data_logchange$Generation=="7",]$sd_logchange^2)
data_logchange[data_logchange$Generation=="7",]$Vector_sample <- data_logchange[data_logchange$Generation=="7",]$Inv_sum /
  sum(data_logchange[data_logchange$Generation=="7",]$Inv_sum,na.rm = TRUE)

data_logchange[data_logchange$Generation=="29",]$Inv_sum <- 1 / (data_logchange[data_logchange$Generation=="29",]$sd_logchange^2)
data_logchange[data_logchange$Generation=="29",]$Vector_sample <- data_logchange[data_logchange$Generation=="29",]$Inv_sum /
  sum(data_logchange[data_logchange$Generation=="29",]$Inv_sum,na.rm = TRUE)

# Problem with NA (if one sd missing): replace by 0 (max value)
data_logchange$Vector_sample[is.na(data_logchange$Vector_sample)] <- 0


##### 
#####  G7
##### 
lm_val_G7 = lm(logchange ~ Treatment + Line + SA + Treatment:Fruit_s, 
               weights = Vector_sample, data = data_logchange[data_logchange$Generation=="7",])

Fratio = anova(lm_val_G7)[3, 3]/anova(lm_val_G7)[4, 3]
pvalue = 1 - pf(Fratio, anova(lm_val_G7)[3, 1], anova(lm_val_G7)[4, 1]) 

pvalue
## [1] 0.9777576
Fratio
## [1] 0.0009158865
anova(lm_val_G7)[3, 1]
## [1] 1
anova(lm_val_G7)[4, 1]
## [1] 3
##### 
#####  G29
##### 
lm_val_G29 = lm(logchange ~ Treatment + Line + SA + Treatment:Fruit_s, 
               weights = Vector_sample, data = data_logchange[data_logchange$Generation=="29",])

Fratio = anova(lm_val_G29)[3, 3]/anova(lm_val_G29)[4, 3]
pvalue = 1 - pf(Fratio, anova(lm_val_G29)[3, 1], anova(lm_val_G29)[4, 1]) 

pvalue
## [1] 0.211062
Fratio
## [1] 2.513403
anova(lm_val_G29)[3, 1]
## [1] 1
anova(lm_val_G29)[4, 1]
## [1] 3
#####  G29
##### 
lm_val_G29 = lm(logchange ~ Treatment  , 
               weights = Vector_sample, data = data_logchange[data_logchange$Generation=="29"&
                                                                data_logchange$SA=="1",])
lm_val_G29_avec = lm(logchange ~ Treatment + Line , 
               weights = Vector_sample, data = data_logchange[data_logchange$Generation=="29"&
                                                                data_logchange$SA=="1",])

anova(lm_val_G29,lm_val_G29_avec)
## Analysis of Variance Table
## 
## Model 1: logchange ~ Treatment
## Model 2: logchange ~ Treatment + Line
##   Res.Df      RSS Df Sum of Sq F Pr(>F)
## 1      8 0.018934                      
## 2      0 0.000000  8  0.018934
summary(lm_val_G29_avec)
## 
## Call:
## lm(formula = logchange ~ Treatment + Line, data = data_logchange[data_logchange$Generation == 
##     "29" & data_logchange$SA == "1", ], weights = Vector_sample)
## 
## Weighted Residuals:
## ALL 11 residuals are 0: no residual degrees of freedom!
## 
## Coefficients: (2 not defined because of singularities)
##                     Estimate Std. Error t value Pr(>|t|)
## (Intercept)           0.1477         NA      NA       NA
## TreatmentCranberry    0.5775         NA      NA       NA
## TreatmentStrawberry   0.2187         NA      NA       NA
## LineCEB               0.1027         NA      NA       NA
## LineCEC               0.2109         NA      NA       NA
## LineCRA              -0.8740         NA      NA       NA
## LineCRB              -0.4279         NA      NA       NA
## LineCRC              -0.6981         NA      NA       NA
## LineCRD              -0.8814         NA      NA       NA
## LineCRE                   NA         NA      NA       NA
## LineFRA              -0.2093         NA      NA       NA
## LineFRB               0.1938         NA      NA       NA
## LineFRC                   NA         NA      NA       NA
## 
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:    NaN 
## F-statistic:   NaN on 10 and 0 DF,  p-value: NA
Fratio = anova(lm_val_G29)[3, 3]/anova(lm_val_G29)[4, 3]
pvalue = 1 - pf(Fratio, anova(lm_val_G29)[3, 1], anova(lm_val_G29)[4, 1]) 

pvalue
## [1] NA
Fratio
## [1] NA
anova(lm_val_G29)[3, 1]
## [1] NA
anova(lm_val_G29)[4, 1]
## [1] NA